*! Source of lmoremata.mlib


*! {smcl}
*! {marker mm_kern}{bf:mm_kern.mata}{asis}
*! version 1.0.3  18aug2020  Ben Jann
version 9.2

local kernels ///
 epanechnikov ///
 epan2 ///
 biweight ///
 triweight ///
 cosine ///
 gaussian ///
 parzen ///
 rectangle ///
 triangle

local klist tokens("`kernels'")

local kname `"(k!="" ? k : "epan2")"'

mata:

string scalar _mm_unabkern(|string scalar k)
{
	return(mm_strexpand(k, `klist', "epan2", 0, k+": invalid kernel"))
}

pointer scalar _mm_findkern(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kern_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kern_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkint(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kint_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kint_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkderiv(|string scalar k)
{
	pointer(real matrix function) scalar K
	K = findexternal("mm_kderiv_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kderiv_"+`kname'+"() not found")
	return(K)
}

pointer scalar _mm_findkdel0(|string scalar k)
{
	pointer(real scalar function) scalar K
	K = findexternal("mm_kdel0_"+`kname'+"()")
	if (K==NULL) _error(3499,"mm_kdel0_"+`kname'+"() not found")
	return(K)
}

real matrix mm_kern(string scalar k, real matrix x)
{
	return((*_mm_findkern(_mm_unabkern(k)))(x))
}

real matrix mm_kint(string scalar k, real scalar r, |real matrix x)
{
	return(mm_callf(mm_callf_setup(_mm_findkint(_mm_unabkern(k)), args()-2, x), r))
}

real matrix mm_kderiv(string scalar k, real matrix x)
{
	return((*_mm_findkderiv(_mm_unabkern(k)))(x))
}

real matrix mm_kdel0(string scalar k)
{
	return((*_mm_findkdel0(_mm_unabkern(k)))())
}

// mm_kern_*: kernel function
// mm_kint_*(1,x): integral of K(z) from -infty to x
// mm_kint_*(2,x): integral of K(z)^2 from -infty to x
// mm_kint_*(2): kernel "roughness"
// mm_kint_*(3,x): integral of z*K(z) from -infty to x
// mm_kint_*(4,x): integral of z^2*K(z) from -infty to x
// mm_kint_*(4): kernel variance
// mm_kvar_*: variance
// mm_kderiv*: derivative of kernel
// mm_kdel0_*: canonical bandwidth

// (scaled) Epanechnikov kernel function
real matrix mm_kern_epanechnikov(real matrix x)
{
	return((.75:-.15*x:^2)/sqrt(5):*(abs(x):<sqrt(5)))
}
real matrix mm_kint_epanechnikov(real scalar r, |real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+(.75*x-.05*x:^3)/sqrt(5)):*(abs(x):<sqrt(5)) + (x:>=sqrt(5)))
	if (r==2 & args()<2) return(3/5^1.5)
	if (r==2)
	 return((.3/sqrt(5):+(.1125*x-.015*x:^3+.0009*x:^5)):*(abs(x):<sqrt(5)) + 3/5^1.5*(x:>=sqrt(5)))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.1875*sqrt(5):+.375/sqrt(5)*x:^2-.0375/sqrt(5)*x:^4):*(abs(x):<sqrt(5)))
	if (r==4 & args()<2) return(1)
	if (r==4)
	 return((.5:+.25/sqrt(5)*x:^3-.03/sqrt(5)*x:^5):*(abs(x):<sqrt(5)) + (x:>=sqrt(5)))
	else _error(3300)
}
//real scalar mm_kvar_epanechnikov() return(1)
real matrix mm_kderiv_epanechnikov(real matrix x)
{
	return(-0.3/sqrt(5)*x:*(abs(x):<sqrt(5)))
}
real scalar mm_kdel0_epanechnikov() return((3/5^1.5)^.2)


// Epanechnikov kernel function
real matrix mm_kern_epan2(real matrix x)
{
	return((.75:-.75*x:^2):*(abs(x):<1))
}
real matrix mm_kint_epan2(real scalar r, |real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.75*x-.25*x:^3):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(.6)
	if (r==2)
	 return((.3:+.5625*x-.375*x:^3+.1125*x:^5):*(abs(x):<1) + .6*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.1875:+.375*x:^2-.1875*x:^4):*(abs(x):<1))
	if (r==4 & args()<2) return(.2) // 1/5
	if (r==4)
	 return((.1:+.25*x:^3-.15*x:^5):*(abs(x):<1) + .2*(x:>=1))
	else _error(3300)
}
//real scalar mm_kvar_epan2() return(.2)
real matrix mm_kderiv_epan2(real matrix x)
{
	return((-1.5*x):*(abs(x):<1))
}
real scalar mm_kdel0_epan2() return(15^.2)


// Biweight kernel function
real matrix mm_kern_biweight(real matrix x)
{
	return(.9375*(1:-x:^2):^2:*(abs(x):<1))
}
real matrix mm_kint_biweight(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.9375*(x-2/3*x:^3+.2*x:^5)):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(5/7)
	if (r==2)
	 return((5/14:+225/256*x-75/64*x:^3+135/128*x:^5-225/448*x:^7+25/256*x:^9):*(abs(x):<1) + 5/7*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.15625:+.46875*x:^2-.46875*x:^4 +.15625*x:^6):*(abs(x):<1))
	if (r==4 & args()<2) return(1/7)
	if (r==4)
	 return((1/14:+.3125*x:^3-.375*x:^5+15/112*x:^7):*(abs(x):<1) + (x:>=1)/7)
	else _error(3300)
}
//real scalar mm_kvar_biweight() return(1/7)
real matrix mm_kderiv_biweight(real matrix x)
{
	return(3.75*(x:^3:-x):*(abs(x):<1))
}
real scalar mm_kdel0_biweight() return(35^.2)


// Triweight kernel function
real matrix mm_kern_triweight(real matrix x)
{
	return(1.09375*(1:-x:^2):^3:*(abs(x):<1))
}
real matrix mm_kint_triweight(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+1.09375*(x-x:^3+.6*x:^5-1/7*x:^7)):*(abs(x):<1) + (x:>=1))
	if (r==2 & args()<2) return(350/429)
	if (r==2)
	 return((175/429:+1225/1024*(x-2*x:^3+3*x:^5-20/7*x:^7+5/3*x:^9-6/11*x:^11+1/13*x:^13)):*(abs(x):<1) + 350/429*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-35/256:+35/64*x:^2-105/128*x:^4+35/64*x:^6-35/256*x:^8):*(abs(x):<1))
	if (r==4 & args()<2) return(1/9)
	if (r==4)
	 return((1/18:+35/96*x:^3-21/32*x:^5+15/32*x:^7-35/288*x:^9):*(abs(x):<1) + (x:>=1)/9)
	else _error(3300)
}
//real scalar mm_kvar_triweight() return(1/9)
real matrix mm_kderiv_triweight(real matrix x)
{
	return(-6.5625*x:*(1:-x:^2):^2:*(abs(x):<1))   // 210/32 = 6.5625
}
real scalar mm_kdel0_triweight() return((9450/143)^.2)


// Cosine trace
real matrix mm_kern_cosine(real matrix x)
{
	return((1:+cos(2*pi()*x)):*(abs(x):<.5))
}
real matrix mm_kint_cosine(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+x+sin(2*pi()*x)/(2*pi())):*(abs(x):<.5) + (x:>=.5))
	if (r==2 & args()<2) return(1.5)
	if (r==2)
	 return((.75:+1.5*x+sin(2*pi()*x)/pi()+cos(2*pi()*x):*sin(2*pi()*x)/(4*pi())):*(abs(x):<.5) + 1.5*(x:>=.5))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.125+.25/pi()^2:+cos(2*pi()*x)/(4*pi()^2)+x:*sin(2*pi()*x)/(2*pi())+.5*x:^2):*(abs(x):<.5))
	if (r==4 & args()<2) return(.5*(1/6-1/pi()^2))
	if (r==4)
	 return((1/24-.25/pi()^2:-sin(2*pi()*x)/(4*pi()^3)+x:*cos(2*pi()*x)/(2*pi()^2)+x:^2:*sin(2*pi()*x)/(2*pi())+x:^3/3):*(abs(x):<.5) + .5*(1/6-1/pi()^2)*(x:>=.5))
	else _error(3300)
}
//real scalar mm_kvar_cosine() return(.5*(1/6-1/pi()^2))
real matrix mm_kderiv_cosine(real matrix x)
{
	return(-2*pi()*sin(2*pi()*x):*(abs(x):<.5))
}
real scalar mm_kdel0_cosine() return((3/(.5*(1/6-1/pi()^2)^2))^.2)


// Gaussian kernel function
real matrix mm_kern_gaussian(real matrix x)
{
	return(normalden(x)) // alternatively: exp(-0.5*((x)^2))/sqrt(2*pi())
}
real matrix mm_kint_gaussian(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)            return(normal(x))
	if (r==2 & args()<2) return(.5/sqrt(pi())) // = normalden(0,sqrt(2))
	if (r==2)            return(.5/sqrt(pi()) * normal(sqrt(2)*x))
	if (r==3 & args()<2) return(0)
	if (r==3)            return(-1/sqrt(2*pi())*exp(-.5*x:^2))
	if (r==4 & args()<2) return(1)
	if (r==4)            return(-1/sqrt(2*pi())*x:*exp(-.5*x:^2) + normal(x))
	else _error(3300)
}
//real scalar mm_kvar_gaussian() return(1)
real matrix mm_kderiv_gaussian(real matrix x)
{
	return(-x:*normalden(x))
}
real scalar mm_kdel0_gaussian() return((1/(4*pi()))^.1)


// Parzen kernel function
real matrix mm_kern_parzen(real matrix x)
{
	return(((4/3:-8*x:^2+8*abs(x):^3):*(abs(x):<=.5)) +
	 ((8*(1:-abs(x)):^3/3):*(abs(x):>1/2:&abs(x):<=1)))
}
real matrix mm_kint_parzen(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((2/3:+8/3*(x+1.5*x:^2+x:^3+.25*x:^4)):*(x:>=-1:&x:<-.5) +
	  (.5:+4/3*x-8/3*x:^3-2*x:^4):*(x:>=-.5:&x:<0) +
	  (.5:+4/3*x-8/3*x:^3+2*x:^4):*(x:>=0:&x:<=.5) +
	  (1/3:+8/3*(x-1.5*x:^2+x:^3-.25*x:^4)):*(x:>.5:&x:<=1) + (x:>1))
	if (r==2 & args()<2) return(302/315)
	if (r==2)
	 return((64/63:+64/9*(x+3*x:^2+5*x:^3+5*x:^4+3*x:^5+x:^6+x:^7/7)):*(x:>=-1:&x:<-.5) +
	  (151/315:+16/9*(x-4*x:^3-3*x:^4+36/5*x:^5+12*x:^6+36/7*x:^7)):*(x:>=-.5:&x:<0) +
	  (151/315:+16/9*(x-4*x:^3+3*x:^4+36/5*x:^5-12*x:^6+36/7*x:^7)):*(x:>=0:&x:<=.5) +
	  (-2/35:+64/9*(x-3*x:^2+5*x:^3-5*x:^4+3*x:^5-x:^6+x:^7/7)):*(x:>.5:&x:<=1) + 302/315*(x:>1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-2/15:+4/3*x:^2+8/3*x:^3+2*x:^4+8/15*x:^5):*(x:>=-1:&x:<-.5) +
	  (-7/60:+2/3*x:^2-2*x:^4-1.6*x:^5):*(x:>=-.5:&x:<0) +
	  (-7/60:+2/3*x:^2-2*x:^4+1.6*x:^5):*(x:>=0:&x:<=.5) +
	  (-2/15:+4/3*x:^2-8/3*x:^3+2*x:^4-8/15*x:^5):*(x:>.5:&x:<=1))
	if (r==4 & args()<2) return(1/12)
	if (r==4)
	 return((2/45:+8/9*x:^3+2*x:^4+1.6*x:^5+4/9*x:^6):*(x:>=-1:&x:<-.5) +
	  (1/24:+4/9*x:^3-1.6*x:^5-4/3*x:^6):*(x:>=-.5:&x:<0) +
	  (1/24:+4/9*x:^3-1.6*x:^5+4/3*x:^6):*(x:>=0:&x:<=.5) +
	  (7/180:+8/9*x:^3-2*x:^4+1.6*x:^5-4/9*x:^6):*(x:>.5:&x:<=1) + (x:>1)/12)
	else _error(3300)
}
//real scalar mm_kvar_parzen() return(1/12)
real matrix mm_kderiv_parzen(real matrix x)
{
	return(((-16*x+24*sign(x):*abs(x):^2):*(abs(x):<=.5)) +
	 ((-8*sign(x):*(1:-abs(x)):^2):*(abs(x):>1/2:&abs(x):<=1)))
}
real scalar mm_kdel0_parzen() return((4832/35)^.2)


// Rectangle kernel function
real matrix mm_kern_rectangle(real matrix x)
{
	return(.5*(round(abs(x),1e-8):<1))
} // 1e-8 same as in kdensity.ado, version 2.3.6 26jun2000 (Stata 7)
real matrix mm_kint_rectangle(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+.5*x):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1))
	if (r==2 & args()<2) return(.5)
	if (r==2)
	 return((.25:+.25*x):*(round(abs(x),1e-8):<1) + .5*(round(x,1e-8):>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-.25:+.25*x:^2):*(round(abs(x),1e-8):<1))
	if (r==4 & args()<2) return(1/3)
	if (r==4)
	 return((1/6:+x:^3/6):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1)/3)
	else _error(3300)
}
//real scalar mm_kvar_rectangle() return(1/3)
real matrix mm_kderiv_rectangle(real matrix x)
{
	return(J(rows(x), cols(x), 0)) // actually, the derivative is infty at -1 and -infty at 1
}
real scalar mm_kdel0_rectangle() return((9/2)^.2)


// Triangle kernel function
real matrix mm_kern_triangle(real matrix x)
{
	return((1:-abs(x)):*(abs(x):<1))
}
real matrix mm_kint_triangle(| real scalar r, real matrix x)
{
	if (r==1 & args()<2) return(1)
	if (r==1)
	 return((.5:+x+.5*x:^2):*(x:>-1:&x:<0) +
	  (.5:+x-.5*x:^2):*(x:>=0:&x:<1) + (x:>=1))
	if (r==2 & args()<2) return(2/3)
	if (r==2)
	 return((x+x:^2+(1:+x:^3)/3):*(x:>-1:&x:<0) +
	  (x-x:^2+(1:+x:^3)/3):*(x:>=0:&x:<1) + 2/3*(x:>=1))
	if (r==3 & args()<2) return(0)
	if (r==3)
	 return((-1/6:+.5*x:^2+x:^3/3):*(x:>-1:&x:<0) +
	  (-1/6:+.5*x:^2-x:^3/3):*(x:>=0:&x:<1))
	if (r==4 & args()<2) return(1/6)
	if (r==4)
	 return((1/12:+x:^3/3+.25*x:^4):*(x:>-1:&x:<0) +
	  (1/12:+x:^3/3-.25*x:^4):*(x:>=0:&x:<1) + (x:>=1)/6)
	else _error(3300)
}
//real scalar mm_kvar_triangle() return(1/6)
real matrix mm_kderiv_triangle(real matrix x)
{
	return(-sign(x):*(abs(x):<1))
}
real scalar mm_kdel0_triangle() return(24^.2)

end


*! {smcl}
*! {marker mm_quantile}{bf:mm_quantile.mata}{asis}
*! version 2.0.6  14jan2022  Ben Jann
version 9.2
mata:

real matrix mm_quantile(real matrix X, | real colvector w,
    real matrix P, real scalar d, real scalar fw, real scalar wd)
{
    real colvector o
    pointer (real matrix) scalar XX, ww
    
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) d = 2
    if (args()<5) fw = 0
    if (!anyof(0::11, d)) {
        display("{err}{it:def} must be an integer in [0,11]")
        _error(3300)
    }
    if (missing(X) | missing(w)) _error(3351)
    if (any(w:<0)) {
        display("{err}{it:w} must not be negative")
        _error(3300)
    }
    // handle weights
    XX = &X; ww = &w
    if (rows(w)!=1) {
        if (rows(w)!=rows(X)) _error(3200)
        // drop zero frequency observations
        if (anyof(w,0)) {
            o = select(1::rows(w), w)
            if (cols(o)==0) o = J(0,1,.) // select() may return 0x0
            XX = &(X[o,])
            ww = &(w[o])
        }
        // weights can be ignored if constant and d < 3
        if (d < 3 & rows(*ww)) { // (*ww may be void)
            if (allof(*ww, (*ww)[1])) ww = &1
        }
        // weights can be ignored if w = 1 for all observations
        else if (allof(*ww,1)) ww = &1
    }
    else if (w==0) {
        XX = &(J(0,cols(X),.))
        ww = &(J(0,1,.))
    }
    // compute quantiles
    if (cols(X)==1 & cols(P)!=1 & rows(P)==1)
        return(_mm_quantile_sort(*XX, *ww, P', d, fw, wd)')
    if (cols(X)!=1 & cols(P)!=1 & cols(X)!=cols(P)) _error(3200)
    return(_mm_quantile_sort(*XX, *ww, P, d, fw, wd))
}

real matrix _mm_quantile_sort(real matrix X, real colvector w,
    real matrix P, real scalar d, real scalar fw, real scalar wd)
{
    real scalar    i, c, c1, c2
    real colvector p, sX, sw, pP, sP
    real matrix    Q

    c1 = cols(X); c2 = cols(P)
    c = max((c1,c2))
    Q = J(rows(P), c, .)
    if (c1==c2) {
        if (w==1) {
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sort(X[,i],1), editmissing(P[,i],1), d, wd)
            }
            return(Q)
        }
        if (rows(w)==1) {
            for (i=c; i; i--) {
                pP = order(P[,i],1); sP = P[pP,i]
                Q[pP,i] = _mm_quantile_w(sort(X[,i],1), w, sP, d, fw, wd)
            }
            return(Q)
        }
        for (i=c; i; i--) {
            p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]
            pP = order(P[,i],1); sP = P[pP,i]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw, wd)
        }
        return(Q)
    }
    if (c1==1) {
        if (w==1) {
            sX = sort(X,1)
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sX, editmissing(P[,i],1), d, wd)
            }
            return(Q)
        }
        if (rows(w)==1) {
            sX = sort(X,1)
            for (i=c; i; i--) {
                pP = order(P[,i],1); sP = P[pP,i]
                Q[pP,i] = _mm_quantile_w(sX, w, sP, d, fw, wd)
            }
            return(Q)
        }
        p = order((X,w),(1,2)); sX = X[p]; sw = w[p]
        for (i=c; i; i--) {
            pP = order(P[,i],1); sP = P[pP,i]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw, wd)
        }
        return(Q)
    }
    if (c2==1) {
        if (w==1) {
            sP = editmissing(P,1)
            for (i=c; i; i--) {
                Q[,i] = _mm_quantile_d(sort(X[,i],1), sP, d, wd)
            }
            return(Q)
        }
        if (rows(w)==1) {
            pP = order(P,1); sP = P[pP]
            for (i=c; i; i--) {
                Q[pP,i] = _mm_quantile_w(sort(X[,i],1), w, sP, d, fw, wd)
            }
            return(Q)
        }
        pP = order(P,1); sP = P[pP]
        for (i=c; i; i--) {
            p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]
            Q[pP,i] = _mm_quantile_w(sX, sw, sP, d, fw, wd)
        }
        return(Q)
    }
    _error(3200)
}

real matrix _mm_quantile(real colvector X, | real colvector w, 
    real matrix P, real scalar d, real scalar fw, real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-negative and non-missing
    real scalar    i
    real colvector o
    real matrix    Q
    pointer (real matrix) scalar XX, ww
    
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) d = 2
    if (args()<5) fw = 0
    // handle weights
    XX = &X; ww = &w
    if (rows(w)!=1) {
        if (rows(w)!=rows(X)) _error(3200)
        // drop zero frequency observations
        if (anyof(w,0)) {
            o = select(1::rows(w), w)
            if (cols(o)==0) o = J(0,1,.) // select() may return 0x0
            XX = &(X[o,])
            ww = &(w[o])
        }
        // weights can be ignored if constant and d < 3
        if (d<3 & rows(*ww)) { // (*ww may be void)
            if (allof(*ww, (*ww)[1])) ww = &1
        }
        // weights can be ignored if w = 1 for all observations
        else if (allof(*ww,1)) ww = &1
    }
    else if (w==0) {
        XX = &(J(0,1,.))
        ww = &(J(0,1,.))
    }
    // compute weighted quantiles: requires sorted p
    if (*ww!=1) {  // not rows(*ww)!=1 !
        if (rows(P)==1) {
            o = order(P',1)
            Q = o // just to dimension q
            Q[o] = _mm_quantile_w(*XX, *ww, P[o]', d, fw, wd)
            return(Q')
        }
        Q = J(rows(P), cols(P), .)
        for (i=cols(P); i; i--) {
            o = order(P[,i],1)
            Q[o,i] = _mm_quantile_w(*XX, *ww, P[o,i], d, fw, wd)
        }
        return(Q)
    }
    // compute unweighted quantiles
    if (rows(P)==1) return(_mm_quantile_d(*XX, editmissing(P',1), d, wd)')
    Q = J(rows(P), cols(P), .)
    for (i=cols(P); i; i--) {
        Q[,i] = _mm_quantile_d(*XX, editmissing(P[,i],1), d, wd)
    }
    return(Q)
}

real colvector _mm_quantile_d(real colvector X, real colvector p, 
    real scalar d, real scalar wd)
{   // X assumed sorted and non-missing
    // p assumed nonmissing
    real scalar     n, eps
    real colvector  pn, j, j1, h

    n = rows(X)
    if ((rows(p)*n)==0) return(J(rows(p), 1, .))   // no obs or rows(p)==0
    if (n==1)           return(J(rows(p), 1, X))   // only one obs
    if (d==10) return(_mm_quantile_d_hd(X, p, wd)) // Harrell-Davis
    if (d==11) return(_mm_quantile_11(X, 1, p))    // mid-quantile by Ma et al.
    if      (d==0) pn = p * n
    else if (d==1) pn = p * n
    else if (d==2) pn = p * n
    else if (d==3) pn = p * n :- .5
                                               // pn = a + p*(n + 1 - a - b)
    else if (d==4) pn = p * n                  //      a = 0, b = 1
    else if (d==5) pn = p * n :+ .5            //      a = b = 0
    else if (d==6) pn = p * (n + 1)            //      a = b = 0 
    else if (d==7) pn = 1   :+ p * (n - 1)     //      a = b = 1
    else if (d==8) pn = 1/3 :+ p * (n + 1/3)   //      a = b = 1/3
    else if (d==9) pn = 3/8 :+ p * (n + 1/4)   //      a = b = 3/8
    else {
        display("{err}{it:def} must be an integer in [0,11]")
        _error(3300)
    }
    if (d==0) return(X[mm_clip(floor(pn):+1, 1, n)])
    if (d==1) return(X[mm_clip(ceil(pn), 1, n)])
    if (d<=3) {
        j = floor(pn)
        if (d==2) h = ((pn:>j) :+ 1) / 2
        else      h = (pn:>j) :| mod(j,2)
    }
    else {
        eps = 4 * epsilon(1) // handle rounding error as in R's quantile()
        j = floor(pn :+ eps)
        h = pn - j
        h = h :* (abs(h):>=eps)
    }
    j1 = mm_clip(j:+1, 1, n)
    _mm_clip(j, 1, n)
    return((1:-h) :* X[j] :+ h :* X[j1])
}

real colvector _mm_quantile_w(real colvector X, real colvector w, real colvector p, 
    real scalar d, real scalar fw, real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-missing and *strictly* positive
    // p assumed sorted
    real scalar n
    real matrix W
    
    n = rows(X)
    if (n==0)       return(J(rows(p), 1, .))
    if (n==1)       return(J(rows(p), 1, X))
    if (rows(w)==0) return(J(rows(p), 1, .))
    // mid-quantile by Ma et al. (2011)
    if (d==11) return(__mm_quantile_11(X, w, p))
    // definition 3 with fw==0
    if (fw==0 & d==3) {
        if (rows(w)==1) W = (1::n) * w
        else {
            W = mm_colrunsum(w, 1, 1) // use quad precision
            if (W[n]>=.) W = J(n,1,.)
        }
        return(_mm_quantile_w_3b(X, W, p))
    }
    // other definitions
    W = _mm_ecdf2(X, w, 0, 1) // W = (uniq X, runningsum(w))
    if (fw==0 & d>3) {
        // definitions 4-10: rescale weights to effective sample size
        if (rows(w)!=1) W[,2] = W[,2] * (W[rows(W),2] / quadsum(w:^2))
    }
    if (d==0)  return(_mm_quantile_w_0(W[,1], W[,2], p))
    if (d==1)  return(_mm_quantile_w_1(W[,1], W[,2], p))
    if (d==2)  return(_mm_quantile_w_2(W[,1], W[,2], p))
    if (d==3)  return(_mm_quantile_w_3(W[,1], W[,2], p))
    if (d==4)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==5)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==6)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==7)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==8)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==9)  return(_mm_quantile_w_d(W[,1], W[,2], p, d))
    if (d==10) return(_mm_quantile_w_hd(W[,1], W[,2], p, wd))
    display("{err}{it:def} must be an integer in [0,11]")
    _error(3300)
}

real colvector _mm_quantile_w_0(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi
    real colvector P, q
    
    j = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<=pi) break
        }
        q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_1(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi
    real colvector P, q
    
    j = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_2(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, r, pi
    real colvector P, q
    
    j = r = rows(W) 
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        if (W[j]==pi) {
            if (j==r) q[i] = x[j]
            else      q[i] = (x[j] + x[j+1])/2
        }
        else q[i] = x[j]
    }
    return(q)
}

real colvector _mm_quantile_w_3(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi, lo, up, d0, d1
    real colvector P, q
    
    j = rows(W)
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        // find j such that W[j-1] < pi <= W[j]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        // reached bottom
        if (j==1) { 
            q[|1\i|] = J(i,1, x[1])
            break
        }
        // if pi is larger than or equal to the next integer after W[j-1],
        // then use the upper x-value; decide between lower and upper x-value
        // only if pi is between W[j-1] and next integer
        lo = W[j-1]
        up = floor(lo) + 1
        if (pi>=up) q[i] = x[j]
        else {
            // set upper boundary to minimum between W[j] and next integer
            // after W[j-1]
            up = min( (W[j], up) )
            // obtain distances between pi and boundaries
            d0 = pi - lo
            d1 = up - pi
            // case 1: pi closer to upper boundary
            if (d1<d0) q[i] = x[j]
            // case 2: pi closer to lower boundary
            else if (d0<d1) q[i] = x[j-1]
            // case 3: equal distance
            else {
                // use x-value from lower boundary if lower boundary is integer 
                // and even; else use x-value from upper boundary; this is an
                // arbitrary decision
                if (!mod(lo,2)) q[i] = x[j-1]
                else            q[i] = x[j]
            }
        }
    }
    return(q)
}

real colvector _mm_quantile_w_3b(real colvector x, real colvector W, real colvector p)
{
    real scalar    i, j, pi, d0, d1
    real colvector P, q
    
    j = rows(W)
    P = p * W[j]
    i = rows(p)
    q = J(i,1,.)
    for (; i; i--) {
        pi = P[i]
        for (; j>1; j--) {
            if (W[j-1]<pi) break
        }
        // reached bottom
        if (j==1) { 
            q[|1\i|] = J(i,1, x[1])
            break
        }
        // obtain distances
        d0 = pi - W[j-1]
        d1 = W[j] - pi
        // case 1: pi closer to upper boundary
        if (d1<d0) q[i] = x[j]
        // case 2: pi closer to lower boundary
        else if (d0<d1) q[i] = x[j-1]
        // case 3: equal distance
        else {
            // use x-value from upper boundary if index of upper boundary is
            // even; else use x-value from lower boundary
            if (!mod(j,2)) q[i] = x[j]
            else           q[i] = x[j-1]
        }
    }
    return(q)
}

real colvector _mm_quantile_w_d(real colvector x, real colvector W, 
    real colvector p, real scalar d) // modifies W
{
    real scalar    N, n, i, r, ll, ul, a, b, hmax
    real colvector h, q, ww
    
    N = rows(x)
    n = W[N]
    if      (d==4) h = p * n
    else if (d==5) h = p * n :+ .5
    else if (d==6) h = p * (n + 1)
    else if (d==7) h = 1   :+ p * (n - 1)
    else if (d==8) h = 1/3 :+ p * (n + 1/3)
    else if (d==9) h = 3/8 :+ p * (n + 1/4)
    r = rows(h)
    if (r) {
        if ((h[1]-1)<=W[1]) {
            // duplicate min if lower limit of first quantile is below data
            x = x[1] \ x 
            W = (h[1]-2) \ W
            N = N + 1
        }
        if (h[r]>W[N]) {
            // duplicate max if upper limit of last quantile is above data
            // (using max(h) instead of h[r] because p may include missing)
            hmax = max(h)
            if (hmax<. & hmax>W[N]) {
                x = x \ x[N]
                W = W \ hmax
                N = N + 1
            }
        }
    }
    q = J(r,1,.)
    a = 2
    for (i=1; i<=r; i++) {
        ul = h[i]; ll = ul - 1
        // exit once lower limit is above data (this also takes care of p>=.)
        if (ll>W[N]) {
            q[|i\.|] = J(r-i+1, 1, x[N])
            break
        }
        // lower index
        while (1) {
            if (W[a]>=ll) break
            if (a>=N) break
            a++
        }
        // upper upper index
        b = a
        if (b<N) {
            while (1) {
                if (W[b]>=ul) break
                if (b>=N) break
                b++
            }
        }
        // compute quantile
        ww = W[|a-1 \ b|] :- ll
        ww = mm_diff(rowmax((J(rows(ww),1,0), rowmin((ww, J(rows(ww),1,1))))))
        q[i] = quadsum(ww :* x[|a \ b|])
    }
    return(q)
}

real colvector _mm_quantile_d_hd(real colvector X, real colvector p,
    real scalar wd)
{
    real scalar    n, i
    real colvector q, W
    
    i = rows(p)
    q = J(i,1,.)
    n = rows(X)
    W = 0 \ (1::n)/n
    for (; i; i--) q[i] = __mm_quantile_hd(X, W, p[i], n, wd)
    return(q)
}

real colvector _mm_quantile_w_hd(real colvector X, real colvector W,
    real colvector p, real scalar wd)  // modifies W
{
    real scalar    n, i
    real colvector q
    
    i = rows(p)
    q = J(i,1,.)
    n = W[rows(W)] 
    W = 0 \ W/n
    for (; i; i--) q[i] = __mm_quantile_hd(X, W, max((0,min((1,p[i])))), n, wd)
    return(q)
}

real scalar __mm_quantile_hd(real colvector X, real colvector W, real scalar p,
    real scalar n, real scalar wd0)
{   // note: domain of a and b in ibeta() is 1e-10 to 1e+17; see help ibeta()
    real scalar    a, b, wd, l, r, il, ir
    real colvector w
    
    a = (n + 1) * p
    b = (n + 1) * (1 - p)
    if      (a<1e-10) return(X[1])
    else if (b<1e-10) return(X[rows(X)])
    if (a>1e+17 | b>1e+17) {
        display("{err}sample size too large; cannot evaluate {bf:ibeta()}")
        _error(3498)
    }
    // trimmed estimator (Akinshin 2021)
    if (wd0<1) {
        wd = wd0<=0 ? 1 / sqrt(n) : wd0
        if      (a<=1 & b<=1) {; l=.   ; r=. ; } // can only happen if n=1
        else if (a<=1 & b>1 ) {; l=0   ; r=wd; } // left border
        else if (a>1  & b<=1) {; l=1-wd; r=1 ; } // right border
        else {; l = __mm_quantile_hd_l(a, b, wd); r = l+wd; }
        if (l<.) { // (else use untrimmed estimator)
            // find lower index
            if (l<=0) il = 1
            else      mm_hunt(W, l, il = floor(l*rows(X)))
            // find upper index
            if (r>=1) ir = rows(W)
            else      {; mm_hunt(W, r, ir = il); ir++; }
            // compute weights
            w = W[|il\ir|]; w[1] = l; w[rows(w)] = r
            w = (ibeta(a, b, w) :- ibeta(a,b,l)) / mm_diff(ibeta(a, b, l\r))
            return(quadsum(mm_diff(w) :* X[|il\ir-1|]))
        }
    }
    // untrimmed estimator
    w = mm_diff(ibeta(a, b, W))
    return(quadsum(w :* X)) // (faster than quadcross(w, X))
}

real scalar __mm_quantile_hd_l(real scalar a, real scalar b, real scalar wd)
{
    real scalar m, lo, up, l, rc
    
    m = (a - 1) / (a + b - 2)
    lo = max((0, m-wd))
    up = min((m, 1-wd))
    rc = mm_root(l=., &___mm_quantile_hd_l(), lo, up, 0, 1000, a, b, wd)
    if (rc) _error(3498, "could not locate highest density interval")
    return(l)
}

real scalar ___mm_quantile_hd_l(real scalar x, real scalar a, real scalar b,
    real scalar wd) return(betaden(a, b, x) - betaden(a, b, x+wd))

real colvector _mm_quantile_11(real colvector X, real colvector w,
    real colvector p)
{   // mid-quantile by Ma et al. (2011)
    // X assumed sorted and non-missing
    // w assumed non-missing and *strictly* positive
    // p not assumed sorted
    real colvector o
    real colvector q
    
    o = order(p,1)
    q = o   // just to dimension q
    q[o] = __mm_quantile_11(X, w, p[o])
    return(q)
}

real colvector __mm_quantile_11(real colvector X, real colvector w,
    real colvector p)
{   // mid-quantile by Ma et al. (2011)
    // X assumed sorted and non-missing
    // w assumed non-missing and *strictly* positive
    // p assumed sorted
    real scalar    i, pi, j, xj, pj
    real colvector q
    real matrix    m
    
    i = rows(p)
    q = J(i,1,.)
    m = _mm_ecdf2(X, w, 1) // mid cdf at unique values: (x, mid cdf)
    j = rows(m)
    // handle p > max(cdf)
    xj = m[j,1]; pj = m[j,2]
    for (; i; i--) {
        if (pj>=p[i]) break
        q[i] = xj
    }
    // handle rest
    for (; i; i--) {
        pi = p[i]
        while (j) {
            pj = m[j,2]
            if (pj<pi) break
            j--
        }
        if (!j) { // reached bottom (p <= min(cdf))
            q[|1\i|] = J(i, 1, m[1,1])
            break
        }
        // note: j must be lower than rows(m) at this point
        xj = m[j,1]
        q[i] = xj + (pi-pj)/(m[j+1,2]-pj) * (m[j+1,1] - xj)
    }
    return(q)
}

end


*! {smcl}
*! {marker mm_median}{bf:mm_median.mata}{asis}
*! version 1.1.3  Ben Jann  14jan2022
version 9.2
mata:

real rowvector mm_median(real matrix X, | real colvector w, real scalar def, 
    real scalar fw, real scalar wd)
{
    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    return(mm_quantile(X, w, .5, def, fw, wd))
}

real scalar _mm_median(real colvector X, | real colvector w, real scalar def, 
    real scalar fw, real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-missing and non-negative
    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    return(_mm_quantile(X, w, .5, def, fw, wd))
}

end


*! {smcl}
*! {marker mm_iqrange}{bf:mm_iqrange.mata}{asis}
*! version 1.1.1  Ben Jann  14jan2022
version 9.2
mata:

real matrix mm_iqrange(real matrix X, | real colvector w, real scalar def, 
    real scalar fw, real scalar wd)
{
    real matrix q

    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    q = mm_quantile(X, w, (.25 \ .75), def, fw, wd)
    return(q[2,] - q[1,])
}

real scalar _mm_iqrange(real colvector X, | real colvector w, real scalar def, 
    real scalar fw, real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-missing and non-negative
    real matrix q

    if (args()<2) w = 1
    if (args()<3) def = 2
    if (args()<4) fw = 0
    q = _mm_quantile(X, w, (.25 \ .75), def, fw, wd)
    return(q[2] - q[1])
}

end


*! {smcl}
*! {marker mm_hdq}{bf:mm_hdq.mata}{asis}
*! version 1.0.1  Ben Jann  14jan2022
version 9.2
mata:

real matrix mm_hdq(real matrix X, | real colvector w, real matrix P,
    real scalar fw, real scalar wd)
{
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) fw = 0
    return(mm_quantile(X, w, P, 10, fw, wd))
}

real matrix _mm_hdq(real colvector X, | real colvector w, real matrix P,
    real scalar fw, real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-negative and non-missing
    if (args()<2) w = 1
    if (args()<3) P = (0, .25, .50, .75, 1)'
    if (args()<4) fw = 0
    return(_mm_quantile(X, w, P, 10, fw, wd))
}

real rowvector mm_hdmed(real matrix X, | real colvector w, real scalar fw,
    real scalar wd)
{
    if (args()<2) w = 1
    if (args()<3) fw = 0
    return(mm_quantile(X, w, .5, 10, fw, wd))
}

real scalar _mm_hdmed(real colvector X, | real colvector w, real scalar fw,
    real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-negative and non-missing
    if (args()<2) w = 1
    if (args()<3) fw = 0
    return(_mm_quantile(X, w, .5, 10, fw, wd))
}

real matrix mm_hdiqr(real matrix X, | real colvector w, real scalar fw,
    real scalar wd)
{
    real matrix q

    if (args()<2) w = 1
    if (args()<3) fw = 0
    q = mm_quantile(X, w, (.25 \ .75), 10, fw, wd)
    return(q[2,] - q[1,])
}

real scalar _mm_hdiqr(real colvector X, | real colvector w, real scalar fw,
    real scalar wd)
{   // X assumed sorted and non-missing
    // w assumed non-negative and non-missing
    real matrix q

    if (args()<2) w = 1
    if (args()<3) fw = 0
    q = _mm_quantile(X, w, (.25 \ .75), 10, fw, wd)
    return(q[2] - q[1])
}

end


*! {smcl}
*! {marker mm_ecdf}{bf:mm_ecdf.mata}{asis}
*! version 1.1.0  10jul2020  Ben Jann
version 9.2
mata:

real matrix mm_ecdf(real matrix X, | real colvector w, real scalar mid, 
    real scalar nonorm, real scalar brk)
{
    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (args()<5) brk = 0
    return(mm_ranks(X, w, (brk!=0 ? 0 : 3), mid, (nonorm!=0 ? 0 : 1)))
}

real colvector _mm_ecdf(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm, real scalar brk)
{
    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (args()<5) brk = 0
    return(_mm_ranks(X, w, (brk!=0 ? 0 : 3), mid, (nonorm!=0 ? 0 : 1)))
}

real matrix mm_ecdf2(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm)
{
    real colvector p

    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0
    if (rows(X)==0) return(J(0,2,.))
    if (rows(w)==1) {
        p = order(X,1)
        return(_mm_ecdf2(X[p], w, mid, nonorm))
    }
    if (rows(w)!=rows(X)) _error(3200)
    p = order((X,w),(1,2))
    return(_mm_ecdf2(X[p], w[p], mid, nonorm))
}

real matrix _mm_ecdf2(real colvector X, | real colvector w, real scalar mid, 
    real scalar nonorm)
{
    real scalar    n
    real colvector p, cdf

    if (args()<2) w = 1
    if (args()<3) mid = 0
    if (args()<4) nonorm = 0

    // compute running sum
    n = rows(X)
    if (rows(w)==1) cdf = (1::n) * w
    else {
        if (rows(w)!=n) _error(3200)
        // treat missing values as missing and use quad precision
        cdf = mm_colrunsum(w, 1, 1)
        if (cdf[n]>=.) cdf = J(n,1,.)
    }

    // remove ties (select last obs in each group)
    p = select(1::n, _mm_unique_tag(X, 1))
    cdf = cdf[p]
    
    // normalize
    if (!nonorm) cdf = cdf / cdf[rows(cdf)]
    
    // midrank adjustment
    if (mid) cdf = cdf :- mm_diff(0\cdf)/2
    
    // return result
    return(X[p], cdf)
}

end


*! {smcl}
*! {marker mm_relrank}{bf:mm_relrank.mata}{asis}
*! version 1.1.0  13jul2020  Ben Jann
version 9.2
mata:

real matrix mm_relrank(
    real matrix X,
    real colvector w,
    real matrix Y,
    | real scalar mid,
      real scalar nonorm, 
      real scalar brk,
      real colvector w2)
{
    if (args()<4) mid = 0
    if (args()<5) nonorm = 0
    if (args()<6) brk = 0
    if (args()<7) w2 = 1
    if (rows(w)!=1 & rows(w)!=rows(X)) _error(3200)
    if (rows(w2)!=1 & rows(w2)!=rows(Y)) _error(3200)
    if (cols(X)==1 & cols(Y)!=1 & rows(Y)==1 & rows(w2)==1)
        return(_mm_relrank_sort(X, w, Y', mid, nonorm, brk, w2)')
    return(_mm_relrank_sort(X, w, Y, mid, nonorm, brk, w2))
}

real matrix _mm_relrank_sort(
    real matrix X,
    real colvector w,
    real matrix Y,
    real scalar mid,
    real scalar nonorm, 
    real scalar brk,
    real colvector w2)
{
    real scalar    i, c, c1, c2
    real colvector p, sX, sw, p2, sY, sw2
    real matrix    R

    c1 = cols(X); c2 = cols(Y)
    c = max((c1,c2))
    R = J(rows(Y), c, .)
    if (c1==c2) {
        for (i=c; i; i--) {
            if (rows(w)==1) {; p = order(X[,i],1); sX = X[p,i]; sw = w; }
            else {; p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]; }
            if (rows(w2)==1) {; p2 = order(Y[,i],1); sY = Y[p2,i]; sw2 = w2; }
            else {; p2 = order((Y[,i],w2),(1,2)); sY = Y[p2,i]; sw2 = w2[p2]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    if (c1==1) {
        if (rows(w)==1) {; p = order(X,1); sX = X[p]; sw = w; }
        else {; p = order((X,w),(1,2)); sX = X[p]; sw = w[p]; }
        for (i=c; i; i--) {
            if (rows(w2)==1) {; p2 = order(Y[,i],1); sY = Y[p2,i]; sw2 = w2; }
            else {; p2 = order((Y[,i],w2),(1,2)); sY = Y[p2,i]; sw2 = w2[p2]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    if (c2==1) {
        if (rows(w2)==1) {; p2 = order(Y,1); sY = Y[p2]; sw2 = w2; }
        else {; p2 = order((Y,w2),(1,2)); sY = Y[p2]; sw2 = w2[p2]; }
        for (i=c; i; i--) {
            if (rows(w)==1) {; p = order(X[,i],1); sX = X[p,i]; sw = w; }
            else {; p = order((X[,i],w),(1,2)); sX = X[p,i]; sw = w[p]; }
            R[p2,i] = _mm_relrank(sX, sw, sY, mid, nonorm, brk, sw2)
        }
        return(R)
    }
    _error(3200)
}

real colvector _mm_relrank(
    real colvector X,
    real colvector w,
    real colvector Y,
    | real scalar mid,
      real scalar nonorm, 
      real scalar brk,
      real colvector w2)
{
    real matrix cdf
    
    if (args()<4) mid = 0
    if (args()<5) nonorm = 0
    if (args()<6) brk = 0
    if (args()<7) w2 = 1
    if (rows(X)==0 | rows(Y)==0) return(J(rows(Y),1,.))
    
    // obtain CDF at unique values of X
    cdf = _mm_ecdf2(X, w, 0, nonorm)
    
    // compute relative ranks
    if (brk==0) {
        if (mid==0) return(_mm_relrank_1(Y, cdf[,1], cdf[,2]))
                    return(_mm_relrank_2(Y, cdf[,1], cdf[,2]))
    }
    if (rows(w2)==1) {
        if (mid==0) return(_mm_relrank_3(Y, cdf[,1], cdf[,2]))
                    return(_mm_relrank_4(Y, cdf[,1], cdf[,2]))
    }
    if (mid==0)     return(_mm_relrank_5(Y, w2, cdf[,1], cdf[,2]))
                    return(_mm_relrank_6(Y, w2, cdf[,1], cdf[,2]))
}

real colvector _mm_relrank_1(real colvector x, real colvector y, real colvector cdf)
{    // case 1: brk = 0, mid = 0
    real scalar    i, j, xi
    real colvector r
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) r[i] = cdf[j]
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_2(real colvector x, real colvector y, real colvector cdf)
{   // case 2: brk = 0, mid = 1
    real scalar    i, j, xi
    real colvector r, step
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) r[i] = cdf[j] - step[j]/2
            else          r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_3(real colvector x, real colvector y, real colvector cdf)
{
    // case 3: brk = 1, mid = 0, no weights
    real scalar    i, j, k, xi
    real colvector r, step

    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            r[i] = cdf[j]
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) continue // no ties
                r[|k\i-1|] = cdf[j] :- step[j] :* (i:-(k::i-1)) / (i-k+1)
                i = k
            }
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_4(real colvector x, real colvector y, real colvector cdf)
{
    // case 4: brk = 1, mid = 1, no weights
    real scalar    i, j, k, xi
    real colvector r, step
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) {
                    r[i] = cdf[j] - step[j] * 0.5
                    continue
                }
                r[|k\i|] = cdf[j] :- step[j] :* ((i+.5):-(k::i)) / (i-k+1)
                i = k
            }
            else r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_5(real colvector x, real colvector w, 
    real colvector y, real colvector cdf)
{
    // case 5: brk = 1, mid = 0, weighted
    real scalar    i, j, k, xi, W
    real colvector r, step, ww
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            r[i] = cdf[j]
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) continue // no ties
                ww = mm_colrunsum(w[|k\i|])
                W  = ww[rows(ww)]
                if (W==0) {
                    r[|k\i-1|] = cdf[j] :- step[j] :* (i:-(k::i-1)) / (i-k+1)
                }
                else {
                    ww = ww[|1\rows(ww)-1|]
                    r[|k\i-1|] = cdf[j] :- step[j] :* (W:-ww) / W
                }
                i = k
            }
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}
real colvector _mm_relrank_6(real colvector x, real colvector w, 
    real colvector y, real colvector cdf)
{
    // case 6: brk = 1, mid = 1, weighted
    real scalar    i, j, k, xi, W
    real colvector r, step, ww
    
    i = rows(x)
    r = J(i, 1, 0)
    j = rows(y)
    step = mm_diff(0\cdf)
    for (; i; i--) {
        xi = x[i]
        for (; j; j--) {
            if (y[j]<=xi) break
        }
        if (j) {
            if (y[j]==xi) {
                for (k=i-1; k; k--) { // find ties in x
                    if (x[k]<xi) break
                }
                if ((++k)==i) {
                    r[i] = cdf[j] - step[j] * 0.5
                    continue // no ties
                }
                ww = mm_colrunsum(w[|k\i|])
                W  = ww[rows(ww)]
                if (W==0) {
                    // if all observations have zero weight
                    r[|k\i|] = cdf[j] :- step[j] :* ((i+.5):-(k::i)) / (i-k+1)
                }
                else {
                    ww = ww - 0.5 * w[|k\i|]
                    r[|k\i|] = cdf[j] :- step[j] :* (W:-ww) / W
                }
                i = k
            }
            else r[i] = cdf[j]
        }
        else break // x[i] is smaller than min(y)
    }
    return(r)
}

end


*! {smcl}
*! {marker mm_ranks}{bf:mm_ranks.mata}{asis}
*! version 1.1.1  30dec2021  Ben Jann
version 9.2
mata:

real matrix mm_ranks(real matrix X, | real colvector w,
 real scalar method, real scalar mid, real scalar norm)
{
    real scalar c, r, i
    real matrix result

    if (args()<2) w = 1
    if (args()<3) method = 0
    if (args()<4) mid = 0
    if (args()<5) norm = 0
    c = cols(X)
    r = rows(X)
    if (rows(w)!=1 & rows(w)!=r) _error(3200)
    if (c<1) return(J(r,0,.))
    if (c==1) return(_mm_ranks_sort(X, w, method, mid, norm))
    result = J(r, c, .)
    for (i=1; i<=c; i++) {
        result[,i] = _mm_ranks_sort(X[,i], w, method, mid, norm)
    }
    return(result)
}

real colvector _mm_ranks_sort(real colvector x, real colvector w,
 real scalar method, real scalar mid, real scalar norm)
{
    real scalar    I
    real colvector p, r

    I = rows(x)
    if (I==0) return(J(0,1,.))
    if      (method==4 & rows(w)!=1) p = order((x,w),(1,2))
    else if (method==0 | method==4)  p = order(x,1)
    else if (rows(w)!=1) p = order((x,w),(1,2)) // include w for stable results
    else p = order(x,1)
    r = __mm_ranks(x[p], (rows(w)!=1 ? w[p] : w), method, mid, norm)
    r[p] = r
    return(r)
}

real colvector _mm_ranks(real colvector x, | real colvector w,
 real scalar method, real scalar mid, real scalar norm) // sorted input assumed
{
    if (args()<2) w = 1
    if (args()<3) method = 0
    if (args()<4) mid = 0
    if (args()<5) norm = 0
    if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
    return(__mm_ranks(x, w, method, mid, norm))
}

real colvector __mm_ranks(real colvector x, real colvector w,
 real scalar method, real scalar mid, real scalar norm) // sorted input assumed
{
    real scalar    i, I, i0, i1
    real colvector ranks, R
    pointer scalar wR
    pragma unset   R

    // compute running sum
    I = rows(x)
    if (I==0) return(J(0,1,.))
    if (rows(w)!=1) {
        // treat missing values as missing and use quad precision
        ranks = mm_colrunsum(w, 1, 1)
        if (ranks[I]>=.) return(J(I,1,.))
    }
    else ranks = (1::I) * w
    
    // no special treatment of ties; use fast computations 
    if (method==0 | method==4) {
        if (mid & norm) return((ranks :- w/2) / ranks[I])
        if (norm)       return(ranks / ranks[I])
        if (mid)        return(ranks :- w/2)
        return(ranks)
    }
    
    // normalize
    if (norm) ranks = ranks / ranks[I]
    
    // handle ties
    if (method==3) {            // use lowest rank
        for(i=I-1; i>=1; i--) {
            if (x[i]==x[i+1]) ranks[i] = ranks[i+1]
        }
    }
    else if (method==1) {       // use highest rank
        for(i=2; i<=I; i++) {
            if (x[i]==x[i-1]) ranks[i] = ranks[i-1]
        }
    }
    else if (method==2) {       // use average rank
        i0 = i1 = 1
        if (rows(w)==1) wR = &1
        else wR = &R
        for(i=2; i<=I; i++) {
            if (x[i0]==x[i]) {
                i1++
                if (i<I) continue
            }
            if (i0<i1) {
                R = (i0 \ i1)
                ranks[|R|] = J(i1-i0+1,1,1) :* mean(ranks[|R|], w[|*wR|])
            }
            i0 = i1 = i
        }
    }
    else _error(3498,"ties = " + strofreal(method) + " not allowed")
    
    // apply midpoint adjustment (at this point ties have same rank)
    if (mid) {
        i0 = ranks[1]
        i1 = i0/2
        ranks[1] = i1
        for (i=2; i<=I; i++) {
            if (x[i]!=x[i-1]) {
                i1 = (i0 + ranks[i])/2
                i0 = ranks[i]
            }
            ranks[i] = i1
        }
    }
    
    // done
    return(ranks)
}

end


*! {smcl}
*! {marker mm_ddens}{bf:mm_ddens.mata}{asis}
*! version 1.0.1  12aug2020  Ben Jann
version 9.2
mata:

real matrix mm_ddens(real colvector X, | real colvector w, real vector minmax0,
    real scalar n0, real scalar h0, real scalar qui) // h0 will be replaced
{
    real scalar    n, N, h, lo, up
    real rowvector minmax
    real colvector AT, W, a
    
    // defaults
    if (args()<2) w  = 1
    if (n0>=.) n0 = 2^14
    if (h0<=0) _error(3300)
    if (args()<6) qui = 0
    
    // check input
    if (missing(X) | missing(w)) _error(3351)
    if (any(w:<0)) {
        display("{err}{it:w} must not be negative")
        _error(3300)
    }
    if (n0<2) _error(3300)
    
    // evaluation grid
    n = 2^ceil(ln(n0)/ln(2)) // round up to next power of 2
    if (length(minmax0)>=1) lo = minmax0[1]
    if (length(minmax0)>=2) up = minmax0[2]
    minmax = minmax(X)
    if (lo>=.) lo = minmax[1] - (minmax[2]-minmax[1])/10 // min - 10% of range
    if (up>=.) up = minmax[2] + (minmax[2]-minmax[1])/10 // min + 10% of range
    if (lo>minmax[1]) {
        display("{err}data smaller than lower bound not allowed")
        _error(3300)
    }
    if (up<minmax[2]) {
        display("{err}data larger than upper bound not allowed")
        _error(3300)
    }
    // bin data
    N  = mm_nobs(X, w)
    AT = rangen(lo, up, n)
    if (mm_issorted(X)) W = _mm_exactbin(X, w, rangen(lo, up, n+1)) / N
    else                W = mm_fastexactbin(X, w, rangen(lo, up, n+1)) / N
        // need to use exact binning because linear binning would introduce
        // some (non-vanishing) bias at the boundaries (doubling the first and
        // last grid count does not seem to help); a consequence of exact 
        // binning is that the density estimate will be slightly shifted/stretched
        // to the left; this error can be substantial if the grid size is small,
        // but it vanishes with increasing grid size
    // obtain discrete cosine transform of binned data
    a = Re( (1 \ 2 * exp(1i * (1::n-1) * pi() / (2*n)))
         :* fft(W[mm_seq(1,n-1,2)] \ W[mm_seq(n,2,2)]) )
    // compute bandwidth (unless provided by user)
    if (h0<.) h = (h0 / (up-lo))^2
    else {
        h = _mm_ddens_h(AT, W, N, (1::n-1):^2, (a[2::n]/2):^2, qui)
        h0 = sqrt(h) * (up-lo) // return optimal bandwidth
    }
    // smooth discrete cosine transform using bandwidth and back-transform
    a = a :* exp(-(0::n-1):^2 * pi()^2 * h/2)
    a = Re(invfft((n * exp(-1i * (0::n-1) * pi() / (2*n))) :* a)) / (up-lo)
    a[mm_seq(1, n, 2) \ mm_seq(2, n, 2)] = a[1::n/2 \ n::n/2+1] // reorder
    // return density estimate and grid
    return((a, AT))
}

real scalar _mm_ddens_h(real colvector AT, real colvector W, real scalar N, 
    real colvector I, real colvector a2, real scalar qui)
{
    real scalar s, hmin, h_os, ax, bx, rc, h
    
    // compute oversmoothed bandwidth
    s = _mm_iqrange(AT, W) / 1.349
    if (s<=0) s = sqrt(variance(AT, W*N))
    else      s = min((s, sqrt(variance(AT, W*N))))
    h_os = (s * (243/(35*N))^.2 * mm_kdel0_gaussian() / (AT[rows(AT)]-AT[1]))^2
    // minimum bandwidth given grid size
    hmin = (.5/(rows(AT)-1) * mm_kdel0_gaussian()/mm_kdel0_rectangle())^2
    // find optimal h
    bx = h_os
    ax = max((hmin, bx/2))
    rc = mm_root(h=., &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
    if (rc==2) { // move down
        bx = ax
        while (1) {
            if (bx<=hmin) break // cannot go below hmin
            ax = max((hmin, bx/2))
            rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
            if (rc==2) bx = ax // continue moving down
            else break // this also stops if there is a change in 
                      //  direction (rc==3) (in this case: h = current bx)
        }
    }
    else if (rc==3) { // move up
        ax = bx; bx = ax*1.5
        while (1) {
            rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
            if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
            else break // this also stops if there is a change in 
                       //  direction (rc==2) (in this case: h = current ax)
        }
    }
    if (h<=hmin) { // solution is smaller than hmin
        ax = h_os; bx = ax*1.5
        rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
        if (rc==2) h = . // moving up does not help
        else if (rc==3) { // continue moving up
            ax = bx; bx = ax*1.5
            while (1) {
                rc = mm_root(h, &_mm_ddens_fixpnt(), ax, bx, 0, 100, N, I, a2)
                if (rc==3) {; ax = bx; bx = ax*1.5; } // continue moving up
                else break // this also stops if there is a change in 
                           //  direction (rc==2) (in this case: h = current ax)
            }
        }
    }
    // return result
    if (h>=.) {
        if (qui==0) display("{txt}(bandwidth estimation failed; using oversmoothed rule)")
        return(h_os)
    }
    return(h)
}

real scalar _mm_ddens_fixpnt(real scalar h, real scalar N, real colvector I, 
    real colvector a2)
{
    real scalar    l, s, K0, c
    real colvector f, t
    
    l = 7
    f = 2 * pi()^(2*l) * sum(I:^l :* a2 :* exp(-I * pi()^2 * h))
    for (s=l-1; s>=2; s--) {
        K0 = mm_prod(mm_seq(1, 2*s-1, 2)) / sqrt(2*pi())
        c  = (1 + (1/2)^(s + 1/2)) / 3
        t  = (2 * c * K0/N/f):^(2/(3 + 2*s))
        f  = 2 * pi()^(2*s) * sum(I:^s :* a2 :* exp(-I * pi()^2 * t))
    }
    return((2 * N * sqrt(pi()) * f)^(-2/5) - h)
}

end



*! {smcl}
*! {marker mm_freq}{bf:mm_freq.mata}{asis}
*! version 1.0.4, Ben Jann, 09oct2020
version 9.1
mata:

real colvector mm_freq(transmorphic matrix x,
 | real colvector w, transmorphic matrix levels)
{
	real colvector p

	if (args()<2) w = 1
	if (args()<3) levels = .
	if (cols(x)==0) return(_mm_freq(x, w, levels))
	if (rows(w)==1) return(_mm_freq(sort(x,1..cols(x)), w, levels))
	p = order(x,1..cols(x))
	return(_mm_freq(x[p,], w[p,], levels))
}

real colvector _mm_freq(transmorphic matrix x,
 | real colvector w, transmorphic matrix levels)
{
	real scalar    i, j, l
	real colvector result

	if (args()<2) w = 1
	if (args()<3) levels = .
	if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
	if (levels==.) levels = _mm_uniqrows(x)
	if (rows(x)==0) return(J(0,1, .))
	l = rows(levels)
	result = J(l,1,0)
	j = 1
	for (i=1; i<=rows(x); i++) {
		for (;j<=l;j++) {
			if (x[i,]<=levels[j,]) break
		}
		if (j>l) break
		if (x[i,]!=levels[j,]) continue
		result[j] = result[j] + (rows(w)!=1 ? w[i] : w)
	}
	return(result)
}

real colvector mm_freq2(transmorphic matrix x,
 | real colvector w)
{
    real colvector p

    if (args()<2) w = 1
    if (cols(x)==0) return(_mm_freq2(x, w))
    p = order(x,1..cols(x))
    if (rows(w)==1) return(_mm_freq2(x[p,],w)[invorder(p)])
    return(_mm_freq2(x[p,],w[p,])[invorder(p)])
}

real colvector _mm_freq2(transmorphic matrix x,
 | real colvector w)
{
    real scalar    i, j
    real colvector result

    if (args()<2) w = 1
    if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200)
    if (rows(x)==0) return(J(0, 1, .))
    result = J(rows(x),1,.)
    j = 1
    for (i=2; i<=rows(x); i++) {
        if (x[i,]!=x[i-1,]) {
            result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ?
                (i-j)*w : sum(w[|j \ i-1|])))
            j = i
        }
    }
    result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ?
        (i-j)*w : sum(w[|j \ i-1|])))
    return(result)
}

end


*! {smcl}
*! {marker mm_histogram}{bf:mm_histogram.mata}{asis}
*! version 1.0.3, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_histogram(
 real colvector x,     // data
 | real colvector w,  // weights
   real colvector g,  // grid (interval borders) (default: 10 intervals)
   real scalar dir)   // 0 right inclusive (default), else left inclusive
{
	real matrix H
	real scalar i

	if (args()<2) w = 1
	if (args()<3) g = rangen(min(x), max(x), 11)
	if (args()<4) dir = 0

	H = J(rows(g)-1, 3, 0)
	for (i=1; i<=rows(H); i++) {
		H[i,2] = g[i+1]-g[i]     // width of bins
		H[i,1] = g[i] + H[i,2]/2 // bin midpoints
	}
	H[.,3] = mm_exactbin(x, w, g, dir) :/ (mm_nobs(x,w) * H[.,2])
	return(H)
}
end


*! {smcl}
*! {marker mm_collapse}{bf:mm_collapse.mata}{asis}
*! version 1.0.2  23apr2021  Ben Jann
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real matrix mm_collapse(
    real matrix    X,              // data
    real colvector w,              // weights (or 1)
    real colvector ID,             // subgroup IDs
    |  pointer(function) scalar f, // function; default: mean()
      `opts')                      // arguments to pass through to f()
{
    real colvector p
    pointer scalar Xs, IDs, ws

    if ((rows(X)!=rows(ID)) | (rows(w)!=1 & rows(X)!=rows(w))) _error(3200)
    p     = order(ID, 1)
    if (isfleeting(X))      Xs  = &(X = X[p,])
    else                    Xs  = &(X[p,])
    if (rows(w)==1)         ws  = &w
    else if (isfleeting(w)) ws  = &(w = w[p])
    else                    ws  = &(w[p])
    if (isfleeting(ID))     IDs = &(ID = ID[p])
    else                    IDs = &(ID[p])

    if (args()==3)  return(_mm_collapse(*Xs, *ws, *IDs))
    if (args()==4)  return(_mm_collapse(*Xs, *ws, *IDs, f))
    if (args()==5)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1))
    if (args()==6)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2))
    if (args()==7)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3))
    if (args()==8)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4))
    if (args()==9)  return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5))
    if (args()==10) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6))
    if (args()==11) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7))
    if (args()==12) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8))
    if (args()==13) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8, o9))
                    return(_mm_collapse(*Xs, *ws, *IDs, f, `opts'))
}


real matrix _mm_collapse(
    real matrix    X,               // data
    real colvector w,               // weights (or 1)
    real colvector ID,              // subgroup IDs
    |  pointer(function) scalar f0, // function; default: mean()
      `opts')                       // arguments to pass through to f()
{
    real scalar                 i, j, c, a, b, ww
    real matrix                 info, res, key
    pointer(function) scalar    f
    transmorphic                setup

    if (rows(X)!=rows(ID))     _error(3200)
    ww      = rows(w)!=1
    if (ww & rows(w)!=rows(X)) _error(3200)
    if (rows(X)<1)             return(J(0, cols(X)+1, missingof(X)))

    f     = args()<4 ? &mean() : f0
    setup = mm_callf_setup(f, args()-4, `opts')

    info = _mm_panels(ID)
    key  = J(rows(info), 1, missingof(ID))
    res  = J(rows(info), c=cols(X), missingof(X))
    b = 0
    for (i=1; i<=rows(info); i++) {
        a = b + 1
        b = b + info[i]
        key[i] = ID[a]
        for (j=1; j<=c; j++) {
            res[i, j] =
                mm_callf(setup, X[|a,j \ b,j|], ww ? w[|a \ b|] : w)
        }
    }
    return(key, res)
}

real matrix mm_collapse2(
    real matrix    X,              // data
    real colvector w,              // weights (or 1)
    real colvector ID,             // subgroup IDs
    |  pointer(function) scalar f, // function; default: mean()
      `opts')                      // arguments to pass through to f()
{
    real colvector p
    pointer scalar Xs, IDs, ws
    real matrix    res

    if ((rows(X)!=rows(ID)) | (rows(w)!=1 & rows(X)!=rows(w))) _error(3200)
    p     = order(ID, 1)
    if (isfleeting(X))      Xs  = &(X = X[p,])
    else                    Xs  = &(X[p,])
    if (rows(w)==1)         ws  = &w
    else if (isfleeting(w)) ws  = &(w = w[p])
    else                    ws  = &(w[p])
    if (isfleeting(ID))     IDs = &(ID = ID[p])
    else                    IDs = &(ID[p])

    res = J(rows(X), cols(X), missingof(X))
    if      (args()==3)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs)
    else if (args()==4)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f)
    else if (args()==5)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1)
    else if (args()==6)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2)
    else if (args()==7)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3)
    else if (args()==8)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4)
    else if (args()==9)  res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5)
    else if (args()==10) res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6)
    else if (args()==11) res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7)
    else if (args()==12) res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8)
    else if (args()==13) res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, o1, o2, o3, o4, o5, o6, o7, o8, o9)
    else                 res[p,] = _mm_collapse2(*Xs, *ws, *IDs, f, `opts')
    return(res)
}

real matrix _mm_collapse2(
    real matrix    X,               // data
    real colvector w,               // weights (or 1)
    real colvector ID,              // subgroup IDs
    |  pointer(function) scalar f0, // function; default: mean()
      `opts')                       // arguments to pass through to f()
{
    real scalar                 i, j, c, a, b, ww
    real matrix                 info, res
    pointer(function) scalar    f
    transmorphic                setup

    if (rows(X)!=rows(ID))     _error(3200)
    ww      = rows(w)!=1
    if (ww & rows(w)!=rows(X)) _error(3200)
    if (rows(X)<1)             return(J(0, cols(X), missingof(X)))

    f     = args()<4 ? &mean() : f0
    setup = mm_callf_setup(f, args()-4, `opts')

    info = _mm_panels(ID)
    res  = J(rows(X), c=cols(X), missingof(X))
    b = 0
    for (i=1; i<=rows(info); i++) {
        a = b + 1
        b = b + info[i]
        for (j=1; j<=c; j++) {
            res[|a,j \ b,j|] = J(info[i], 1,
                mm_callf(setup, X[|a,j \ b,j|], ww ? w[|a \ b|] : w))
        }
    }
    return(res)
}

end


*! {smcl}
*! {marker mm_gini}{bf:mm_gini.mata}{asis}
*! version 1.0.3  26mar2008  Ben Jann
version 9.2
mata:

real rowvector mm_gini(real matrix X, | real colvector w)
{
    real rowvector result
    real scalar c, i

    if (args()==1) w = 1
    c = cols(X)
    if (c<1) return(J(1,0,.))
    if (c==1) return(_mm_gini(X, w))
    result = J(1, c, .)
    for (i=1; i<=c; i++) {
        result[1, i] = _mm_gini(X[,i], w)
    }
    return(result)
}

real scalar _mm_gini(real colvector x, real colvector w)
{
    real matrix mv

    if (rows(x)<1) return(.)
    mv = mm_meanvariance0((x, mm_ranks(x, w, 3, 1, 1)), w)
    return(mv[3,1] * 2 / mv[1,1])
}

end


*! {smcl}
*! {marker mm_nobs}{bf:mm_nobs.mata}{asis}
*! version 1.0.0, Ben Jann, 07jun2006
version 9.2
mata:

real scalar mm_nobs(transmorphic matrix x, real colvector w)
{
	return(rows(w)==1 ? rows(x)*w :
	 (rows(x)==rows(w) ? colsum(w) : _error(3200)))
}
end


*! {smcl}
*! {marker mm_greedy}{bf:mm_greedy.mata}{asis}
*! version 1.0.1  17may2019  Ben Jann
version 9.2
local Int   real scalar
local IntC  real colvector
local IntM  real matrix
local RS    real scalar
local RC    real colvector
local RM    real matrix
local T     transmorphic
local TM    transmorphic matrix
local PidC  pointer(`Int') colvector
local PidM  pointer(`Int') matrix
local dist  pointer(function) scalar
mata:

`RM' mm_greedy2(`TM' T, `TM' C, `Int' n, `RS' calip, `dist' f, | `T' fopts)
{
    return(mm_greedy_pairs(mm_greedy(T, C, n, calip, f, fopts)))
}

`RM' mm_greedy_pairs(`IntM' P)
{
    `Int'  i, j, k, n
    `IntC' N
    `RM'   E
    
    N = rownonmissing(P)
    k = sum(N)
    E = J(k, 3, .)
    for (i=rows(P); i; i--) {
        n = N[i]
        for (j=n; j; j--) {
            E[k,] = (i, P[i,j], 1/n)
            k--
        }
    }
    return(E)
}

`IntM' mm_greedy(`TM' T, `TM' C, `Int' n, `RS' calip, `dist' f, | `T' fopts)
{
    `RC'   d
    `PidM' ij
    pragma unset d
    
    // check input
    if (cols(T)!=cols(C)) _error(3200, "T and C must have same number of columns")
    if (cols(T)==0) _error(3200, "T and C must have at least one column")
    // make ID lookup table and compute distances
    if (calip>=.) ij = _mm_greedy_dist(d, T, C, f, fopts)
    else          ij = _mm_greedy_dist_calip(d, T, C, f, fopts, calip)
    // match and return IDs of matched controls
    if (n<=1 | n>=.) return(_mm_greedy_match_1(ij, d, rows(T)))
    return(_mm_greedy_match_n(ij, d, rows(T), n))
}

`PidM' _mm_greedy_dist(`RC' d, `TM' T, `TM' C, `dist' f, `T' fopts)
{
    `Int'  i, a, b, nT, nC, N
    `PidC' Ti, Ci
    `PidM' ij
    
    nT = rows(T); nC = rows(C)   // number of treated and controls
    N = nT * nC                  // number of combinations
    Ti = _mm_greedy_pid(nT)      // make ID pointers for treated
    Ci = _mm_greedy_pid(nC)      // make ID pointers for controls
    d = J(N, 1, .)               // prepare vector of distances
    ij = J(N, 2, NULL)           // prepare ID lookup table
    if (N==0) return(ij)         // T or C is void
    b = 0
    for (i=nT; i; i--) {         // loop over treated => compare to all controls
        a = b + 1                // start index of current batch
        b = b + nC               // end index of current batch
        d[|a\b|] = (*f)(T[i,], C, fopts)   // add distances to distance vector
        ij[|a,1\b,.|] = J(nC,1,Ti[i]), Ci  // add ID pointers to lookup table
    }
    return(ij)
}

`PidM' _mm_greedy_dist_calip(`RC' d, `TM' T, `TM' C, `dist' f, `T' fopts,
    `RS' calip)
{
    `Int'  i, a, b, nT, nC, N, n
    `PidC' Ti, Ci
    `PidM' ij
    `RC'   D
    `IntC' P, p
    
    nT = rows(T); nC = rows(C)    // number of treated and controls
    N = nT * nC                   // number of combinations
    Ti = _mm_greedy_pid(nT)       // make ID pointers for treated
    Ci = _mm_greedy_pid(nC)       // make ID pointers for controls
    d = J(N, 1, .)                // prepare vector of distances
    ij = J(N, 2, NULL)            // prepare ID lookup table
    if (N==0) return(ij)          // T or C is void
    b = 0; P = (1::nC)
    for (i=nT; i; i--) {          // loop over treated => compare to all controls
        D = (*f)(T[i,], C, fopts) // compute distances
        p = select(P, D:<=calip)  // index of distances within caliper
        n = length(p)             // number of distances within caliper
        if (n==0) continue        // no valid distances -> skip case
        a = b + 1                 // start index of current batch
        b = b + n                 // end index of current batch
        if (n==nC) {              // all distances are valid
            d[|a\b|] = D          // add distances to distance vector
            ij[|a,1\b,.|] = J(nC,1,Ti[i]), Ci // add ID pointers to lookup table
            continue
        }
        d[|a\b|] = D[p]          // add distances to distance vector
        ij[|a,1\b,.|] = J(n,1,Ti[i]), Ci[p] // add ID pointers to lookup table
    }
    if (b==N) return(ij)         // all distances valid; all rows filled
    if (b==0) {                  // no valid distances
        d = J(0, 1, .)
        return(J(0, 2, NULL))
    }
    d = d[|1\b|]                 // remove unused rows from d
    return(ij[|1,1\b,.|])        // remove unused rows from ij
}

`PidC' _mm_greedy_pid(`Int' n)
{
    `Int'  i
    `PidC' P
    
    P = J(n, 1, NULL)                 // prepare container for ID pointers
    for (i=n; i; i--) P[i] = &(i*1)   // generate ID pointers
    return(P)
}

`IntC' _mm_greedy_match_1(`PidM' ij, `RC' d, `Int' nT)
{
    `Int'  i, j, ii, jj, N, L
    `IntC' p, M
    
    M = J(nT, 1, .)                        // prepare container for matched IDs
    L = nT                                 // maximum number of matches
    N = rows(d)                            // number of distances
    p = order(d, 1)                        // obtain sort order of distances
    for (i=1; i<=N; i++) {
        j = p[i]                           // get next position
        if ((ii = *(ij[j,1]))==.) continue // treatment case already matched
        if ((jj = *(ij[j,2]))==.) continue // control case already used
        M[ii] = jj                         // copy matched control ID
        *(ij[j,1]) = .                     // delete treatment ID
        *(ij[j,2]) = .                     // delete control ID
        if (!(--L)) break                  // maximum reached; no treatment IDs left
    }
    return(M)
}

`IntM' _mm_greedy_match_n(`PidM' ij, `RC' d, `Int' nT, `Int' n)
{
    `Int'  i, j, ii, jj, N, L, nc
    `IntC' p, NC
    `IntM' M
    
    NC = J(nT, 1, 0)                       // counter for number of matches
    M = J(nT, n, .)                        // container for matched IDs
    L = nT * n                             // maximum number of matches
    N = rows(d)                            // number of distances
    p = order(d, 1)                        // obtain sort order of distances
    for (i=1; i<=N; i++) {
        j = p[i]                           // get next position
        if ((ii = *(ij[j,1]))==.) continue // treatment case already matched
        if ((jj = *(ij[j,2]))==.) continue // control case already used
        nc = NC[ii]; nc++                  // update number of matches
        M[ii, nc] = jj                     // copy matched control ID
        if (nc==n) *(ij[j,1]) = .          // delete treatment ID
        *(ij[j,2]) = .                     // delete control ID
        NC[ii] = nc                        // store number of matches
        if (!(--L)) break                  // maximum reached; no treatment IDs left
    }
    return(M)
}

end




*! {smcl}
*! {marker mm_colvar}{bf:mm_colvar.mata}{asis}
*! version 1.0.1, Ben Jann, 22jun2006
version 9.0
mata:

real rowvector mm_colvar(real matrix X, |real colvector w)
{
	real rowvector result
	real scalar j, k

	if (args()==1) w = 1

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = variance(X[,j], w)
	}
	return(result)
}

real matrix mm_meancolvar(real matrix X, |real colvector w)
{
	real matrix result
	real scalar j, k

	if (args()==1) w = 1

	result = J(2, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[.,j] = meanvariance(X[,j], w)
	}
	return(result)
}

end


*! {smcl}
*! {marker mm_variance0}{bf:mm_variance0.mata}{asis}
*! version 1.0.1, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_variance0(real matrix X, |real colvector w)
{
	real rowvector	CP
	real rowvector	means
	real scalar	n

	if (args()==1) w = 1

	CP = quadcross(w,0, X,1)
	n  = cols(CP)
	means = CP[|1\n-1|] :/ CP[n]
	if (missing(CP) | CP[n]==0) return(J(cols(X), cols(X), .))
	return(crossdev(X,0,means, w, X,0,means) :/ CP[n])
}

real matrix mm_meanvariance0(real matrix X, |real colvector w)
{
	real rowvector	CP
	real rowvector	means
	real scalar	n

	if (args()==1) w = 1

	CP = quadcross(w,0, X,1)
	n  = cols(CP)
	means = CP[|1\n-1|] :/ CP[n]
	if (missing(means)) return(means \ J(cols(X),cols(X),.))
	return(means \ crossdev(X,0,means, w, X,0,means) :/ CP[n])
}

end


*! {smcl}
*! {marker mm_mse}{bf:mm_mse.mata}{asis}
*! version 1.0.0, Ben Jann, 22jun2006
version 9.0
mata:

real matrix mm_mse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector CP

	CP = cross(w,0, X,1)
	if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .))
	return(mm_sse(X, w, mu) :/ CP[cols(CP)])
}

real matrix mm_sse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector CP

	CP = cross(w,0, X,1)
	if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .))
	return(crossdev(X,0,mu, w, X,0,mu))
}

real rowvector mm_colmse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector result
	real scalar j, k

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = mm_mse(X[,j], w, mu[j])
	}
	return(result)
}

real rowvector mm_colsse(real matrix X, real colvector w, real rowvector mu)
{
	real rowvector result
	real scalar j, k

	result = J(1, k=cols(X), .)
	for (j=1; j<=k; j++) {
		result[j] = mm_sse(X[,j], w, mu[j])
	}
	return(result)
}

end


*! {smcl}
*! {marker mm_sample}{bf:mm_sample.mata}{asis}
*! version 1.0.2, Ben Jann, 23oct2020
version 9.1
mata:

real colvector mm_sample(
 real colvector n,
 real matrix strata,
 | real colvector cluster,
   real colvector w,
   real scalar wor,
   real scalar count,
   real scalar fast,
   real scalar alt,
   real scalar nowarn)
{
	real scalar i, b, e, n0, n1, N, cb, ce, ups
	real colvector s, si, ci, wi, nn, R

	if (args()<3) cluster=.
	if (args()<4) w=1
	if (args()<5) wor=0
	if (args()<6) count=0
	if (args()<7) fast=0
	if (args()<8) alt=0
	if (args()<9) nowarn=0
	if (cols(strata)<1) _error(3200)
	ups = (rows(w)!=1)
	if (fast==0) {
		if (rows(strata)>1 | (cluster==. & ups==0)) {
			if (missing(strata)) _error(3351, "'strata' has missing values")
		}
	}
	if (rows(n)<1) _error(3200)

//check w
	if (ups) {
		if (cluster!=.) {
			if (rows(w)!=rows(cluster))
				_error(3200, "rows('w') unequal number of clusters")
		}
		if (fast==0) {
			if (missing(w)) _error(3351, "'w' has missing values")
			if (colsum(w:<0)) _error(3498, "'w' has negative values")
			if (colsum(w)<=0) _error(3498, "sum('w') is zero")
		}
	}

//simple sample, unstratified
	if (rows(strata)==1) {
		N = strata[1,1]
		if (cluster==.) {
			if (N>=.) N = rows(w)
			else if (ups) {
				if (N!=rows(w))
					_error(3200, "rows('w') unequal population size")
			}
			return(_mm_sample(n, (ups ? w : N), ups, wor, count, alt, nowarn))
		}

//cluster sample, unstratified
		if (rows(cluster)<1) _error(3200)
		if (N>=.) N = colsum(cluster)
		else if (fast==0) {
			if (N!=colsum(cluster))
			  _error(3200, "sum of cluster sizes unequal population size")
		}
		s = _mm_sample(n, (ups ? w : rows(cluster)), ups, wor, count, alt, nowarn)
		return(_mm_expandclusters(s, cluster, count, N))
	}

//stratified sample
	if (rows(strata)<1) _error(3200)
	if (cluster!=. & cols(strata)<2) _error(3200)
//-sample sizes
	if (rows(n)==1) {
		if (n>=.) {
			if (cluster==.) nn = strata[.,1]
			else            nn = strata[.,2]
		}
		else nn = J(rows(strata),1,n)
	}
	else {
		if (rows(n)!=rows(strata)) _error(3200)
		nn = n
		for (i=1;i<=rows(nn);i++) {
			if (nn[i]>=.) {
				if (cluster==.) nn[i] = strata[i,1]
				else            nn[i] = strata[i,2]
			}
		}
	}
//-prepare sample vector
	if (count) s = J(colsum(strata[.,1]),1,.)
	else s = J(colsum(nn),1,.)
//-draw samples within strata
	if (count==0) n0 = 1
	e = 0
	ce = 0
	for (i=1;i<=rows(strata);i++) {
		b = e + 1
		e = e + strata[i,1]
//--simple sample
		if (cluster==.) {
			si = _mm_sample(nn[i], (ups ? w[|b \ e|] : strata[i,1]), ups, wor,
				count, alt, nowarn)
		}
//--cluster sample
		else {
			cb = ce + 1
			ce = ce + strata[i,2]
			ci = cluster[|cb \ ce|]
			wi = (ups ? w[|cb \ ce|] : rows(ci))
			if (fast==0) {
				if (strata[i,1]!=colsum(ci))
				  _error(3200, "sum of cluster sizes unequal size of stratum")
			}
			si = _mm_sample(nn[i], wi, ups, wor, count, alt, nowarn)
			if (count) si = _mm_expandclusters(si, ci, 1, strata[i,1])
		}
//--add subsample to sample vector
		if (count) R = (b \ e)
		else {
			n1 = n0 + rows(si) - 1
			R = (n0 \ n1)
			if (cluster==.) si = si :+ (b-1)
			else si = si :+ (cb-1)
			n0 = n1 + 1
		}
		if (R[1]<=R[2]) s[|R|] = si
	}
	if (count==0&cluster!=.) s = _mm_expandclusters(s, cluster)
	return(s)
}

real colvector _mm_sample(real scalar n, real colvector NorW, 
    real scalar ups, real scalar wor, real scalar count, real scalar alt, 
    real scalar nowarn)
{
    if (ups) {
        if (wor) return(mm_upswor(n, NorW, count, nowarn))
        return(mm_upswr(n, NorW, count))
    }
    if (wor) return(mm_srswor(n, NorW, count, alt))
    return(mm_srswr(n, NorW, count))
} 

real colvector _mm_expandclusters(real colvector s0,
 real colvector cluster, | real scalar count, real scalar N)
{
	real scalar i, e, b, n0, n1
	real colvector s, eclust

	if (args()<3) count = 0
	if (count) {
		s = J(N,1,.)
		e = 0
		for (i=1;i<=rows(cluster);i++) {
			b = e + 1
			e = e + cluster[i]
			s[|b \ e|] = J(e-b+1, 1, s0[i])
		}
		return(s)
	}
	s = J(colsum(cluster[s0,]),1,.)
	eclust = mm_colrunsum(cluster)
	n0 = 1
	for (i=1;i<=rows(s0);i++) {
		e = eclust[s0[i]]
		b = e - cluster[s0[i]] + 1
		n1 = n0 + e-b
		s[|n0 \ n1|] = (b::e)
		n0 = n1+1
	}
	return(s)
}

end


*! {smcl}
*! {marker mm_srswr}{bf:mm_srswr.mata}{asis}
*! version 1.0.0, Ben Jann, 27mar2006
version 9.1
mata:

// simple random sample, with replacement
real colvector mm_srswr(real scalar n, real scalar N,
 | real scalar count)
{
	real colvector res, u
	real scalar i, nn

// check args
	if (args()<3) count = 0
	if (N>=.) _error(3351)
	if (N<1) _error(3300)
	if (n>=.) nn = N
	else nn = n

// no sampling
	if (N==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	u = ceil(uniform(nn,1)*N)
	if (count==0) return(u)
	res = J(N,1,0)
	for (i=1;i<=nn;i++) {
		res[u[i]] = res[u[i]] + 1
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_srswor}{bf:mm_srswor.mata}{asis}
*! version 1.0.1, Ben Jann, 23oct2020
version 9.1
mata:

// simple random sample, without replacement
real colvector mm_srswor(real scalar n, real scalar N,
 | real scalar count, real scalar alt)
{
	real colvector res, u
	real scalar i, nn

// check args
	if (args()<3) count = 0
	if (args()<4) alt = 0
	if (N>=.) _error(3351)
	if (N<1) _error(3300)
	if (n<0) _error(3300)
	if (n>=.) nn = N
	else nn = n
	if (N<nn) _error(3300, "n may not be larger than N")

// no sampling
	if (N==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}
	if (nn==0) {
		if (count) return(J(N,1,0))
		return(J(0,1,.))
	}

// draw sample
	if (alt) {
		if (nn^2 < N/2) u = _mm_srswor_a(n, N)
		else            u = _mm_srswor_b(n, N)
	}
	else u = mm_unorder2(N)[|1 \ nn|] //=> stable results even for very large N
	if (count==0) return(u)
	res = J(N,1,0)
	for (i=1;i<=nn;i++) {
		res[u[i]] = 1
	}
	return(res)
}

real colvector _mm_srswor_a(real scalar n, real scalar N)
{   // naive algorithm; fast when n is small compared to N
    real scalar    i, u
    real colvector H
    
    H = J(n, 1, .)
    H[1] = ceil(uniform(1,1)*N)
    for (i=2; i<=n; i++) {
        while (anyof(H[|1 \ i|], u = ceil(uniform(1,1)*N))) continue
        H[i] = u
    }
    return(H)
}

real colvector _mm_srswor_b(real scalar n, real scalar N)
{   // Fisher–Yates shuffle; fast when n is large compared to N
    real scalar    i, j, k
    real colvector I, H
    
    k = N
    I = 1::k
    H = J(n, 1, .)
    for (i=1; i<=n; i++) {
        j = ceil(uniform(1,1)*k)
        H[i] = I[j]
        I[j] = I[k]
        k--
    }
    return(H)
}


end


*! {smcl}
*! {marker mm_upswr}{bf:mm_upswr.mata}{asis}
*! version 1.0.0, Ben Jann, 01apr2006
version 9.1
mata:

// unequal probability sampling, with replacement
real colvector mm_upswr(real scalar n, real colvector w,
 | real scalar count)
{
	real colvector ub, res, u
	real scalar i, j, nn

// check args
	if (args()<3) count = 0
	if (rows(w)<1) _error(3498, "no cases")
	if (n>=.) nn = rows(w)
	else nn = n

// no sampling
	if (rows(w)==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	ub = mm_colrunsum(w)
	ub = ub/ub[rows(ub)]
	u = sort(uniform(nn,1),1)
	j=1
	if (count) {
		res = J(rows(w),1,0)
		for (i=1;i<=nn;i++) {
			while (u[i]>ub[j]) j++
			res[j] = res[j]+1
		}
		return(res)
	}
	res = J(nn,1,0)
	for (i=1;i<=nn;i++) {
		while (u[i]>ub[j]) j++
		res[i] = j
	}
	return(mm_jumble2(res)) //=> stable results even for very large n
}

end


*! {smcl}
*! {marker mm_upswor}{bf:mm_upswor.mata}{asis}
*! version 1.0.1, Ben Jann, 12jun2007
version 9.1
mata:

// unequal probability sampling, without replacement
real colvector mm_upswor(real scalar n, real colvector w,
 | real scalar count, real scalar nowarn)
{
	real colvector p, ub, res
	real scalar i, j, nn, u, z

// check args
	if (args()<4) nowarn = 0
	if (args()<3) count = 0
	if (rows(w)<1) _error(3498, "no cases")
	if (n<0) _error(3300)
	if (n>=.) nn = rows(w)
	else nn = n
	if (nowarn==0) {
		if (rows(w)<nn) _error(3300, "n may not be larger than rows(w)")
	}

// no sampling
	if (rows(w)==1) {
		if (count) return(nn)
		return(J(nn,1,1))
	}

// draw sample
	p = mm_unorder2(rows(w)) //=> stable results even for very large N
	ub = w[p] / colsum(w) * nn
	if (nowarn==0) {
		z = colsum(ub:>1)
		if (z==1) _error(3300, "1 case has w_i*n/sum(w)>1")
		if (z>0) _error(3300, strofreal(z) + " cases have w_i*n/sum(w)>1")
	}
	ub = mm_colrunsum(ub)
	u = uniform(1,1)
	j=1
	if (count) {
		res = J(rows(w),1,0)
		for (i=1;i<=nn;i++) {
			while (u>ub[j]) j++
			res[j] = res[j]+1
			u = u + 1
		}
		return(res[invorder(p)])
	}
	res = J(nn,1,0)
	for (i=1;i<=nn;i++) {
		while (u>ub[j]) j++
		res[i] = p[j]
		u = u + 1
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_jk}{bf:mm_jk.mata}{asis}
*! version 1.0.0, Ben Jann, 04jul2006
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

struct mm_jkstats {
	real colvector reps, ndrop, fpc
	real rowvector stat
	real matrix    rstat
}

// note: reps is total n. of replications
//       ndrop is n. of unsuccessful replications
//       rows(rstat) = reps-ndrop is n. of successful replications

struct mm_jkstats scalar mm_jk(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar nodots,       // 0
   real colvector strata,    // .
   real colvector cluster,   // .
   real colvector subpop,    // .
   real colvector fpc,       // .
   real matrix stat,         // .
   `opts')                   // opts to pass to f
{
	real scalar               i, ii, i0, i1, h, h0, h1, j
	real colvector            C, wjk, wjki
	real matrix               S
	transmorphic              fsetup
	struct mm_jkstats scalar  jk
	pointer scalar            ws

	if (args()<3) w = 1
	if (args()<4) nodots = 0
	if (args()<7) subpop = .
	if (args()<8) fpc = .
	if (args()<9) stat = .

	if (subpop!=. & rows(subpop)>0) ws = &(w :* (subpop:!=0))
	else ws = &w
	if (rows(*ws)==1) ws = &J(rows(X), 1, *ws)

	fsetup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) jk.stat = mm_callf(fsetup, X, *ws)[1,]
	else         jk.stat = stat[1,]

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")

	if (nodots==0) {
		printf("{txt}\nJackknife replications ({res}%g{txt})\n",
		 colsum(S[,2]))
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	jk.rstat = J(colsum(S[,2]), cols(jk.stat), .)
	jk.ndrop = J(rows(S), 1, 0)
	jk.fpc   = J(rows(S), 1, 0)
	h1 = i1 = 0
	i = ii = 1
	for (h=1; h<=rows(S); h++) { // for each stratum
		h0 = h1 + 1
		h1 = h1 + S[h,1]
		if (fpc!=. & rows(fpc)>0) jk.fpc[h] = fpc[h1]
		wjk = *ws
		wjk[| h0 \ h1 |] = wjk[| h0 \ h1 |] * (S[h,2] / (S[h,2] - 1))
		for (j=1; j<=S[h,2]; j++) { // for each PSU
			i0 = i1 + 1
			i1 = i1 + (C!=. ? C[i] : 1)
			wjki = wjk
			wjki[| i0 \ i1 |] = J(i1-i0+1,1,0)
			jk.rstat[ii,] = mm_callf(fsetup, X, wjki)[1,]
			if (missing(jk.rstat[ii,])) {
				jk.ndrop[h] = jk.ndrop[h] + 1
				if (nodots==0) printf("{err}x{txt}")
			}
			else {
				ii++
				if (nodots==0) printf(".")
			}
			if (nodots==0) {
				if (mod(i,50)==0) printf(" %5.0f\n",i)
				displayflush()
			}
			i++
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) jk.rstat = jk.rstat[|1,1 \ ii-1,.|]
		else jk.rstat = J(0, cols(jk.rstat), .)
	}
	jk.reps = S[,2]
	return(jk)
}

real matrix mm_jk_report(
 struct mm_jkstats scalar jk,    // jackknife estimates
 | string vector what,           // b,pseudo,mean,bias,v,se,ci
   real scalar level,            //
   real scalar mse,              // !=0 => use mse formula
   real vector fpc0)             // sampling fraction (for each stratum)
{
	real scalar    theta, pseudo, mean, bias, se, v, ci, l,
	               i, alpha, tval
	real rowvector jkse, jkv, jkmean
	real colvector fpc, nh, mh
	real matrix    res, jkpseudo

	if (args()<2) what  = "se"
	if (args()<3) level = st_numscalar("c(level)")
	if (args()<4) mse   = 0

	theta = pseudo = mean = bias = v = se = ci = l = 0
	for (i=1; i<=length(what); i++) {
		if      (what[i]=="theta" | what[i]=="b") theta = (l = l + 1)
		else if (what[i]=="pseudo") pseudo = (l = l + rows(jk.rstat))
		else if (what[i]=="mean"  ) mean  = (l = l + 1)
		else if (what[i]=="bias"  ) bias  = (l = l + 1)
		else if (what[i]=="v"     ) v     = (l = l + cols(jk.rstat))
		else if (what[i]=="se"    ) se    = (l = l + 1)
		else if (what[i]=="ci"    ) ci    = (l = l + 2)
		else _error(3498, "'" + what[i] + "' not allowed")
	}

	res = J(l, cols(jk.rstat), .)
	nh = jk.reps - jk.ndrop
	fpc = fpc0
	if (cols(fpc)!=1) fpc = fpc'
	if (rows(fpc)==0 | fpc==.) fpc = jk.fpc
	if (rows(fpc)==0)          fpc = J(rows(jk.reps),1,0)
// observed parameter vector
	if (theta) res[theta,] = jk.stat
// pseudo values
	if (pseudo | (mse==0 & (mean | bias))) {
		jkpseudo = jk.rstat +
		 mm_expand(nh, nh, 1, 1):*(jk.stat :- jk.rstat)
		if (pseudo)
		 res[|pseudo-rows(jk.rstat)+1,1 \ pseudo,.|] = jkpseudo
	}
// jackknife mean
	if (mean | bias) {
		if (mse) jkmean = mean(jk.rstat, 1)
		else     jkmean = mean(jkpseudo, 1)
		if (mean) res[mean,] = jkmean
	}
// bias
	if (bias) res[bias,] = jkmean - jk.stat
// variance
	if (v) {
		mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1)
		if (mse) jkv = mm_sse(jk.rstat, mh, jk.stat)
		else     jkv = mm_sse(jk.rstat, mh, mean(jk.rstat,1))
		res[|v-cols(res)+1,1 \ v,.|] = jkv
	}
// standard errors
	if (se | ci) {
		if (v) jkse = sqrt(diagonal(jkv))'
		else {
			mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1)
			if (mse) jkse = sqrt(mm_colsse(jk.rstat, mh, jk.stat))
			else     jkse = sqrt(mm_colsse(jk.rstat, mh, mean(jk.rstat,1)))
		}
		if (se) res[se,] = jkse
	}
// confidence interval
	if (ci) {
		alpha = (100-level)/200
		tval  = invttail(rows(jk.rstat)-rows(jk.reps), alpha)
		res[ci-1,] = jk.stat - tval*jkse
		res[ci,]   = jk.stat + tval*jkse
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_bs}{bf:mm_bs.mata}{asis}
*! version 1.0.0, Ben Jann, 10jul2006
version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

struct mm_bsstats {
	real scalar    reps, ndrop
	real rowvector stat, se
	real matrix    rstat, rse
}

// note: reps is total n. of replications
//       ndrop is n. of unsuccessful replications
//       rows(rstat) = reps-ndrop is n. of successful replications

struct mm_bsstats scalar mm_bs(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar reps,        // 50
   real scalar d,           // 0
   real scalar nodots,      // 0
   real colvector strata,   // .
   real colvector cluster,  // .
   real matrix stat,        // .
   `opts')                  // opts to pass to f
{
	real scalar               i, ii, hasse
	real colvector            C, N, s
	real matrix               S, res
	transmorphic              setup
	struct mm_bsstats scalar  bs

	if (args()<3) w = 1
	if (args()<4) reps = 50
	if (args()<5) nodots = 0
	if (args()<6) d = 0
	if (args()<9) stat = .

	setup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) {
		res = mm_callf(setup, X, w)
		bs.stat = res[1,]
		if (hasse=(rows(res)>1)) bs.se = res[2,]
	}
	else {
		bs.stat = stat[1,]
		if (hasse=(rows(stat)>1)) bs.se = stat[2,]
	}

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")
	if (d==0|d>=.) N = .
	else N = S[,2] :- d

	if (nodots==0) {
		printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps)
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	bs.rstat = J(reps, cols(bs.stat), .)
	if (hasse) bs.rse = J(reps, cols(bs.se), .)
	for (i=ii=1; i<=reps; i++) {
		s  = mm_sample(N, S, C)
		res = mm_callf(setup, X[s,], rows(w)==1 ? w : w[s,])
		if (missing(res)) {
			if (nodots==0) printf("{err}x{txt}")
		}
		else {
			bs.rstat[ii,] = res[1,]
			if (hasse) bs.rse[ii,] = res[2,]
			ii++
			if (nodots==0) printf(".")
		}
		if (nodots==0) {
			if (mod(i,50)==0) printf(" %5.0f\n",i)
			displayflush()
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) {
			bs.rstat = bs.rstat[|1,1 \ ii-1,.|]
			if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|]
		}
		else {
			bs.rstat = J(0, cols(bs.rstat), .)
			if (hasse) bs.rse = J(0, cols(bs.rse), .)
		}
	}
	bs.reps = reps
	bs.ndrop = reps - rows(bs.rstat)
	return(bs)
}

struct mm_bsstats scalar mm_bs2(
 pointer(real rowvector function) scalar f,
 real matrix X,
 | real colvector w,
   real scalar reps,        // 50
   real scalar d,           // 0
   real scalar nodots,      // 0
   real colvector strata,   // .
   real colvector cluster,  // .
   real matrix stat,        // .
   `opts')                  // opts to pass to f
{
	real scalar               i, ii, hasse
	real colvector            C, N, ws
	real matrix               S, res
	transmorphic              setup
	struct mm_bsstats scalar  bs

	if (args()<3) w = 1
	if (args()<4) reps = 50
	if (args()<5) nodots = 0
	if (args()<6) d = 0
	if (args()<9) stat = .

	setup = mm_callf_setup(f, args()-9, `opts')

	if (stat==.) {
		res = mm_callf(setup, X, w)
		bs.stat = res[1,]
		if (hasse=(rows(res)>1)) bs.se = res[2,]
	}
	else {
		bs.stat = stat[1,]
		if (hasse=(rows(stat)>1)) bs.se = stat[2,]
	}

	mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.)
	if (anyof(S, 1)) _error(460, "singleton cluster detected")
	if (d==0|d>=.) N = .
	else N = S[,2] :- d

	if (nodots==0) {
		printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps)
		display("{txt}{hline 4}{c +}{hline 3} 1 " +
		 "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " +
		 "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ")
	}
	bs.rstat = J(reps, cols(bs.stat), .)
	if (hasse) bs.rse = J(reps, cols(bs.se), .)
	for (i=ii=1; i<=reps; i++) {
		ws = mm_sample(N, S, C, 1, 0, 1):*w
		res = mm_callf(setup, X, ws)
		if (missing(res)) {
			if (nodots==0) printf("{err}x{txt}")
		}
		else {
			bs.rstat[ii,] = res[1,]
			if (hasse) bs.rse[ii,] = res[2,]
			ii++
			if (nodots==0) printf(".")
		}
		if (nodots==0) {
			if (mod(i,50)==0) printf(" %5.0f\n",i)
			displayflush()
		}
	}
	if (nodots==0 & mod(i-1,50)) display("")
	if (ii<i) {
		if (ii>1) {
			bs.rstat = bs.rstat[|1,1 \ ii-1,.|]
			if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|]
		}
		else {
			bs.rstat = J(0, cols(bs.rstat), .)
			if (hasse) bs.rse = J(0, cols(bs.rse), .)
		}
	}
	bs.reps = reps
	bs.ndrop = reps - rows(bs.rstat)
	return(bs)
}

real matrix mm_bs_report(
 struct mm_bsstats scalar bs, // boostrap estimates
 | string vector what,        // b,mean,bias,v,se,n,basic,p,bc,bca,t
   real scalar level,         //
   real scalar mse,           // !=0 => use mse formula
   struct mm_jkstats jk)      // jackknife estimates
{
	real scalar    theta, mean, bias, v, se, n, basic, p, bc,
	               bca, t, l, i, alpha, tval
	real rowvector bsv, bsse, bsmean, z0, a
	real matrix    res

	if (args()<2) what  = "se"
	if (args()<3) level = st_numscalar("c(level)")
	if (args()<4) mse   = 0

	theta = mean = bias = v = se = n = basic = p = bc = bca = t = l = 0
	for (i=1; i<=length(what); i++) {
		if      (what[i]=="theta" | what[i]=="b") theta = (l = l + 1)
		else if (what[i]=="mean" ) mean  = (l = l + 1)
		else if (what[i]=="bias" ) bias  = (l = l + 1)
		else if (what[i]=="v"    ) v     = (l = l + cols(bs.rstat))
		else if (what[i]=="se"   ) se    = (l = l + 1)
		else if (what[i]=="n" | what[i]=="ci") n = (l = l + 2)
		else if (what[i]=="basic") basic = (l = l + 2)
		else if (what[i]=="p"    ) p     = (l = l + 2)
		else if (what[i]=="bc"   ) bc    = (l = l + 2)
		else if (what[i]=="bca"  ) bca   = (l = l + 2)
		else if (what[i]=="t"    ) t     = (l = l + 2)
		else _error(3498, "'" + what[i] + "' not allowed")
	}

	res = J(l, cols(bs.rstat), .)
// observed parameter vector
	if (theta) res[theta,] = bs.stat
// mean of bs replicates
	if (mean | bias) {
		bsmean = mean(bs.rstat, 1)
		if (mean) res[mean,] = bsmean
	}
// bias
	if (bias) res[bias,] = bsmean - bs.stat
// variance
	if (v) {
		if (mse) bsv = mm_mse(bs.rstat, 1, bs.stat)
		else     bsv = variance(bs.rstat, 1)
		res[|v-cols(res)+1,1 \ v,.|] = bsv
	}
// standard errors
	if (se | n) {
		if (v) bsse = sqrt(diagonal(bsv))'
		else {
			if (mse) bsse = sqrt(mm_colmse(bs.rstat, 1, bs.stat))
			else     bsse = sqrt(mm_colvar(bs.rstat, 1))
		}
		if (se) res[se,] = bsse
	}
// normal ci
	alpha = (100-level)/200
	if (n) {
		tval = invttail(bs.reps-bs.ndrop-1, alpha)
		res[n-1,] = bs.stat - tval*bsse
		res[n,]   = bs.stat + tval*bsse
	}
// basic ci
	if (basic) res[|basic-1,1 \ basic,.|] = 2*bs.stat :-
	 mm_quantile(bs.rstat, 1, ((1-alpha) \ alpha))
// percentile ci
	if (p) res[|p-1,1 \ p,.|] =
	 mm_quantile(bs.rstat, 1, (alpha \ (1-alpha)))
// bias-corrected ci
	if (bc | bca) z0 = editmissing(invnormal(
	 colsum(bs.rstat:<=bs.stat)/rows(bs.rstat) ), 0)
	if (bc) {
		res[|bc-1,1 \ bc,.|] = mm_quantile(bs.rstat, 1,
		 (normal(2*z0 :+ invnormal(alpha)) \
		 normal(2*z0 :- invnormal(alpha))))
	}
// bias-corrected and accelerated ci
	if (bca) {
		a = colsum((mean(jk.rstat, 1):-jk.rstat):^3) :/
		    ( 6 * colsum((mean(jk.rstat, 1):-jk.rstat):^2):^1.5 )
		res[|bca-1,1 \ bca,.|] = mm_quantile(bs.rstat, 1,
		 (normal(z0 :+ (z0 :+ invnormal(alpha)) :/
		   (1:-a:*(z0 :+ invnormal(alpha)))) \
		 normal(z0 :+ (z0 :- invnormal(alpha)) :/
		   (1:-a:*(z0 :- invnormal(alpha))))))
	}
// percentile-t ci
	if (t) {
		res[|t-1,1 \ t,.|] = bs.stat :- bs.se :* mm_quantile(
		 editmissing((bs.rstat:-bs.stat):/bs.rse, 0), 1,
		 ((1-alpha) \ alpha) )
	}
	return(res)
}

end


*! {smcl}
*! {marker mm_ls}{bf:mm_ls.mata}{asis}
*! version 1.0.1  22mar2021  Ben Jann
version 9.2
mata:

real colvector mm_lsfit(real colvector y, | real matrix X, real colvector w,
    real scalar cons, real scalar qd, real scalar demean)
{
    if (args()<3) w = 1
    return(mm_ls_b(mm_ls(y, X, w, cons, qd, demean)))
}

struct mm_ls_struct {
    pointer(real colvector) scalar y, w
    pointer(real matrix) scalar    X
    real scalar    cons
    real scalar    k    // number of predictors
    real scalar    kadj // number of non-collinear predictors
    real scalar    N    // sum of weights
    real scalar    ymean
    real rowvector means
    real scalar    rss // residual sum of squares
    real scalar    s   // scale (RMSE)
    real scalar    r2  // R-squared
    real matrix    SSinv, XXinv, V
    real colvector omit
    real colvector b
}

struct mm_ls_struct scalar mm_ls(
    real colvector   y,
    | real matrix    X,
      real colvector w,
      real scalar    cons,   // cons=0 excludes the constant
      real scalar    qd,     // qd=0 uses double precision for cross products
      real scalar    demean) // demean=0 does not use demeaning
{
    real scalar     ymean
    real rowvector  means
    real colvector  ybar, p, Xy, beta
    real matrix     Xbar, XX
    struct mm_ls_struct scalar t
    
    // initialize
    if (args()<3) w = 1
    t.y = &y
    t.X = &X
    t.w = &w
    t.cons = (cons!=0)
    t.k = cols(X)
    t.N = t.ymean = t.means = t.rss = t.s = t.r2 = t.XXinv = t.V = .z
    if (!cons & !t.k) { // model has no parameters
        t.b = t.omit = J(0,1,.)
        t.kadj = 0
        t.XXinv = J(0,0,.)
        t.rss = .
        return(t)
    }
    if (rows(y)==0) {
        t.b = t.omit = J(t.k+t.cons,1,.)
        t.XXinv = J(t.k+t.cons,t.k+t.cons,.)
        t.rss = .
        return(t)
    }
    
    // without constant
    if (cons==0) {
        XX = qd ? quadcross(X, w, X) : cross(X, w, X)
        t.XXinv = invsym(XX)
        t.omit = diagonal(t.XXinv):==0
        p = select(1::t.k, !t.omit)
        t.kadj = length(p)
        if (t.kadj==0) {
            t.b = J(t.k,1,0)
            return(t)
        }
        if (t.kadj<t.k) {
            XX = XX[p,p]
            Xy = qd ? quadcross(X[,p], w, y) : cross(X[,p], w, y)
        }
        else Xy = qd ? quadcross(X, w, y) : cross(X, w, y)
        beta = lusolve(XX, Xy)
        t.b = J(t.k, 1, 0)
        t.b[p] = beta
        return(t)
    }
    
    // constant only
    if (t.k==0) {
        t.kadj = t.omit = 0
        t.b = mm_ls_ymean(t)
        return(t)
    }
    
    // without demeaning
    if (!demean) {
        XX = qd ? quadcross(X,1, w, X,1) : cross(X,1, w, X,1)
        t.XXinv = invsym(XX, t.k+1) // do not omit constant
        t.omit = diagonal(t.XXinv)[|1\t.k|]:==0
        p = select(1::t.k, !t.omit)
        t.kadj = length(p)
        if (t.kadj==0) { // all collinear
            t.b = J(t.k,1,0) \ mm_ls_ymean(t) // add constant
            t.omit = t.omit \ 0
            return(t)
        }
        if (t.kadj<t.k) {
            XX = XX[p \ t.k+1, p \ t.k+1]
            Xy = qd ? quadcross(X[,p],1, w, y,0) : cross(X[,p],1, w, y,0)
        }
        else Xy = qd ? quadcross(X,1, w, y,0) : cross(X,1, w, y,0)
        beta = lusolve(XX, Xy)
        t.b = J(t.k+1, 1, 0)
        t.b[p \ t.k+1] = beta
        t.omit = t.omit \ 0 // add constant
        return(t)
    }
    
    // mean deviation method
    ymean = mm_ls_ymean(t)
    means = mm_ls_means(t)
    Xbar = X :- means
    XX = qd ? quadcross(Xbar, w, Xbar) : cross(Xbar, w, Xbar)
    t.SSinv = invsym(XX)
    t.omit = diagonal(t.SSinv):==0
    p = select(1::t.k, !t.omit)
    t.kadj = length(p)
    if (t.kadj==0) { // all collinear
        t.b = J(t.k,1,0) \ ymean // add constant
        t.omit = t.omit \ 0
        return(t)
    }
    if (t.kadj<t.k) {
        Xbar  = Xbar[,p]
        XX    = XX[p,p]
        means = means[p]
    }
    ybar = y :- ymean
    Xy = qd ? quadcross(Xbar, w, ybar) : cross(Xbar, w, ybar)
    beta = lusolve(XX, Xy)
    t.b = J(t.k+1, 1, 0)
    t.b[p]   = beta
    t.b[t.k+1] = (ymean - means*beta) // add constant
    t.omit = t.omit \ 0
    return(t)
}

real colvector mm_ls_b(struct mm_ls_struct scalar t)
{
    return(t.b)
}

real colvector mm_ls_omit(struct mm_ls_struct scalar t)
{
    return(t.omit)
}

real colvector mm_ls_k_omit(struct mm_ls_struct scalar t)
{
    return(t.k - t.kadj)
}

real matrix mm_ls_V(struct mm_ls_struct scalar t)
{
    if (t.V==.z) t.V = mm_ls_XXinv(t)*mm_ls_s(t)^2
    return(t.V)
}

real colvector mm_ls_se(struct mm_ls_struct scalar t)
{
    return(sqrt(diagonal(mm_ls_V(t))))
}

real scalar mm_ls_N(struct mm_ls_struct scalar t)
{
    if (t.N==.z) t.N = rows(*t.w)==1 ? *t.w * rows(*t.y) : quadsum(*t.w)
    return(t.N)
}

real scalar mm_ls_ymean(struct mm_ls_struct scalar t)
{
    if (t.ymean==.z) {
        if (rows(*t.y)==0) t.ymean = .
        else               t.ymean = quadcross(*t.w, *t.y) / mm_ls_N(t)
    }
    return(t.ymean)
}

real rowvector mm_ls_means(struct mm_ls_struct scalar t)
{
    if (t.means==.z) {
        if      (t.k==0)        t.means = J(1, 0, .)
        else if (rows(*t.X)==0) t.means = J(1, t.k, .)
        else                    t.means = quadcross(*t.w, *t.X) / mm_ls_N(t)
    }
    return(t.means)
}

real matrix mm_ls_XXinv(struct mm_ls_struct scalar t)
{
    // note: XXinv will already be filled in case of cons=0 or demean=0
    if (t.XXinv==.z) {
        if (t.k) {
            t.XXinv = t.SSinv \ -(t.means * t.SSinv)
            t.XXinv = (t.XXinv, (t.XXinv[t.k+1,]' \
                       1/mm_ls_N(t) - t.means * t.XXinv[t.k+1,]'))
        }
        else t.XXinv = 1/mm_ls_N(t) // constant-only model
    }
    return(t.XXinv)
}

real colvector mm_ls_xb(struct mm_ls_struct scalar t, | real matrix X)
{
    if (args()==2) {
        if (cols(X)!=t.k) _error(3200)
        if (t.k==0) {
            if (t.cons) return(J(rows(X), 1, t.b))
            return(J(rows(X), 1, .)) // model has no parameters
        }
        if (t.cons) return(X * t.b[|1\t.k|] :+ t.b[t.k+1])
        return(X * t.b)
    }
    if (t.k==0) {
        if (t.cons) return(J(rows(*t.y), 1, t.b))
        return(J(rows(*t.y), 1, .)) // model has no parameters
    }
    if (t.cons) return(*t.X * t.b[|1\t.k|] :+ t.b[t.k+1])
    return(*t.X * t.b)
}

real scalar mm_ls_rss(struct mm_ls_struct scalar t)
{
    if (t.rss==.z) {
        t.rss = quadcross(*t.w, (*t.y - mm_ls_xb(t)):^2)
    }
    return(t.rss)
}

real scalar mm_ls_s(struct mm_ls_struct scalar t)
{
    if (t.s==.z) {
        t.s = sqrt(mm_ls_rss(t) * editmissing(1 / (mm_ls_N(t)-t.kadj-t.cons), 0))
    }
    return(t.s)
}

real scalar mm_ls_r2(struct mm_ls_struct scalar t)
{
    if (t.r2==.z) {
        t.r2 = 1 - mm_ls_rss(t)/quadcross(*t.w, (*t.y :- mm_ls_ymean(t)):^2)
    }
    return(t.r2)
}

end


*! {smcl}
*! {marker mm_areg}{bf:mm_areg.mata}{asis}
*! version 1.0.4  15sep2021  Ben Jann
version 9.2
mata:

// Main functions

real colvector mm_aregfit(real colvector y, real colvector id, | real matrix X,
    real colvector w, real scalar sort, real scalar qd)
{
    if (args()<3) X = .
    if (args()<4) w = 1
    return(mm_areg_b(mm_areg(y, id, X, w, sort, qd)))
}

struct mm_areg_struct {
    pointer(real colvector) scalar y, w
    pointer(real matrix) scalar    X
    real scalar                    qd, k, a, N, ymean, rss, s, r2
    real colvector                 yd, xb, e, u
    real rowvector                 means
    real matrix                    Xd, XXinv, V
    struct mm_ls_struct scalar     ls
    pointer(struct mm_areg_struct_g scalar) scalar g
}

struct mm_areg_struct_g {
    real scalar    sort             // 1 if data requires sorting, 0 else
    real colvector p                // sort index; only filled in if sort==1
    real colvector idx, levels, n
}

struct mm_areg_struct scalar mm_areg(
    real colvector   y,
    real colvector   id,
    | real matrix    X,
      real colvector w,
      real scalar    sort,
      real scalar    qd)
{
    if (args()<3) X = .
    if (args()<4) w = 1
    if (rows(id)!=rows(y)) _error(3200)
    return(_mm_areg(_mm_areg_g(id, sort), y, X, w, qd))
}

struct mm_areg_struct_g scalar _mm_areg_g(real colvector id, real scalar sort)
{
    struct mm_areg_struct_g scalar g
    
    g.sort = (sort!=0)
    if (g.sort) {
        g.p = mm_order(id, 1, 1) // stable sort order
        __mm_areg_g(g, id[g.p])
    }
    else {
        __mm_areg_g(g, id)
    }
    return(g)
}

void __mm_areg_g(struct mm_areg_struct_g scalar g, real colvector id)
{
    real scalar r
    
    r = rows(id)
    if (r==0) {
        g.idx = g.levels = g.n = J(0,1,.)
    }
    else if (r==1) {
        g.idx = 1
        g.levels = id
        g.n = r
    }
    else {
        g.idx = (id :!= (id[|2\.|] \ id[1]))
        g.idx[r] = 1 // should id be constant
        g.idx = select(1::r, g.idx) // index of last obs in each group
        g.levels = id[g.idx]
        g.n = g.idx - (0 \ g.idx)[|1 \ rows(g.idx)|] // group sizes
    }
}

struct mm_areg_struct scalar _mm_areg(
    struct mm_areg_struct_g scalar g,
    real colvector                 y,
    | real matrix                  X,
      real colvector               w,
      real scalar                  qd)
{
    struct mm_areg_struct scalar t
    
    // setup
    if (args()<3) X = .
    if (args()<4) w = 1
    t.qd = (qd!=0)
    t.a = t.xb = t.e = t.u = .z
    t.N = t.rss = t.s = t.r2 = t.XXinv = t.V = t.ymean = t.means = .z
    t.g = &g
    t.y = &y
    if (X==.) t.X = &(J(rows(y), 0, .))
    else      t.X = &X
    t.w = &w
    t.k = cols(*t.X)
    
    // transform data
    _mm_areg_demean(t.yd, t.Xd, y, *t.X, w, g, t.qd)
    
    // estimate
    t.ls = mm_ls(t.yd, t.Xd, w, 0, t.qd)
    
    // return
    return(t)
}

void _mm_areg_demean(real colvector yd, real matrix Xd, 
    real colvector y, real matrix X, real colvector w,
    struct mm_areg_struct_g scalar g, | real scalar qd)
{
    real matrix M
    
    M = (y,X)
    _mm_areg_gmean(g, M, w, qd)
    yd = y - M[,1]
    if (cols(X)) Xd = X - M[,2..cols(M)]
    else         Xd = X
}

void _mm_areg_gmean(struct mm_areg_struct_g scalar g, real matrix X,
    real colvector w, | real scalar qd)
{
    if (g.sort) {
        _collate(X, g.p)
        if (qd) __mm_areg_gmean(X, rows(w)==1 ? w : w[g.p], g.idx)
        else    __mm_areg_gmean_dbl(X, rows(w)==1 ? w : w[g.p], g.idx)
        _collate(X, invorder(g.p))
    }
    else {
        if (qd) __mm_areg_gmean(X, w, g.idx)
        else    __mm_areg_gmean_dbl(X, w, g.idx)
    }
}

void __mm_areg_gmean(real matrix X, real colvector w, real colvector idx)
{
    real scalar    i, n, a, b, k, W
    real colvector ww
    real matrix    ab
    
    if (rows(X)<1) return
    n = rows(idx)
    b = 0
    if (rows(w)==1) {
        for (i=1; i<=n; i++) {
            a = b + 1
            b = idx[i]
            k = b - a
            if (!k) continue // no averaging necessary if less than two obs
            k++
            ab = (a,1 \ b,.)
            X[|ab|] = J(k, 1, quadcolsum(X[|ab|])/k)
        }
        return
    }
    for (i=1; i<=n; i++) {
        a = b + 1
        b = idx[i]
        k = b - a
        if (!k) continue // no averaging necessary if less than two obs
        k++
        ab = (a,1 \ b,.)
        ww = w[|ab|]
        W  = quadsum(ww)
        if (W) X[|ab|] = J(k, 1, quadcross(ww/W, X[|ab|]))
        else   X[|ab|] = J(k, 1, quadcolsum(X[|ab|])/k)
            // using unweighted average if sum of weights is 0
    }
}

void __mm_areg_gmean_dbl(real matrix X, real colvector w, real colvector idx)
{   // double-precision variant of __mm_areg_gmean()
    real scalar    i, n, a, b, k, W
    real colvector ww
    real matrix    ab

    if (rows(X)<1) return
    n = rows(idx)
    b = 0
    if (rows(w)==1) {
        for (i=1; i<=n; i++) {
            a = b + 1
            b = idx[i]
            k = b - a
            if (!k) continue // no averaging necessary if less than two obs
            k++
            ab = (a,1 \ b,.)
            X[|ab|] = J(k, 1, colsum(X[|ab|])/k)
        }
        return
    }
    for (i=1; i<=n; i++) {
        a = b + 1
        b = idx[i]
        k = b - a
        if (!k) continue // no averaging necessary if less than two obs
        k++
        ab = (a,1 \ b,.)
        ww = w[|ab|]
        W  = sum(ww)
        if (W) X[|ab|] = J(k, 1, cross(ww/W, X[|ab|]))
        else   X[|ab|] = J(k, 1, colsum(X[|ab|])/k)
            // using unweighted average if sum of weights is 0
    }
}

// Results

real colvector mm_areg_b(struct mm_areg_struct scalar t)
{
    return(mm_areg_beta(t) \ mm_areg_alpha(t))
}

real colvector mm_areg_beta(struct mm_areg_struct scalar t)
{
    return(mm_ls_b(t.ls))
}

real scalar mm_areg_alpha(struct mm_areg_struct scalar t)
{
    if (t.a==.z) {
        if (t.k) t.a = mm_areg_ymean(t) - mm_ls_xb(t.ls, mm_areg_means(t))
        else     t.a = mm_areg_ymean(t)
    }
    return(t.a)
}

real colvector mm_areg_xb(struct mm_areg_struct scalar t, | real matrix X)
{
    if (args()==2) {
        if (cols(X)!=t.k) _error(3200)
        if (t.k) return(mm_ls_xb(t.ls, X) :+ mm_areg_alpha(t))
        return(J(rows(X), 1, mm_areg_alpha(t)))
    }
    if (t.xb==.z) {
        if (t.k) t.xb = mm_ls_xb(t.ls, *t.X) :+ mm_areg_alpha(t)
        else     t.xb = J(rows(*t.y), 1, mm_areg_alpha(t))
    }
    return(t.xb)
}

real colvector mm_areg_ue(struct mm_areg_struct scalar t)
{
    return(*t.y - mm_areg_xb(t))
}

real colvector mm_areg_xbu(struct mm_areg_struct scalar t)
{
    return(mm_areg_xb(t) + mm_areg_u(t))
}

real colvector mm_areg_u(struct mm_areg_struct scalar t)
{
    if (t.u==.z) {
        t.u = *t.y - mm_areg_xb(t)
        _mm_areg_gmean(*t.g, t.u, *t.w, t.qd)
        // could also computed t.u = *t.y - mm_areg_e(t) - mm_areg_xb(t), but
        // would still have to take care of roundoff error to make sure that
        // t.u is truly constant within group
    }
    return(t.u)
}

real colvector mm_areg_e(struct mm_areg_struct scalar t)
{
    if (t.e==.z) {
        if (t.k) t.e = t.yd - mm_ls_xb(t.ls)
        else     t.e = t.yd
    }
    return(t.e)
}

real scalar mm_areg_s(struct mm_areg_struct scalar t)
{
    if (t.s==.z) {
        t.s = sqrt(mm_areg_rss(t) * 
            editmissing(1 / (mm_areg_N(t) - (t.k - mm_ls_k_omit(t.ls))
                 - rows(t.g->idx)), 0))
    }
    return(t.s)
}

real scalar mm_areg_r2(struct mm_areg_struct scalar t)
{
    if (t.r2==.z) {
        t.r2 = 1 - mm_areg_rss(t) / 
                   quadcross(*t.w, (*t.y :- mm_areg_ymean(t)):^2)
    }
    return(t.r2)
}

real colvector mm_areg_se(struct mm_areg_struct scalar t)
{
    return(sqrt(diagonal(mm_areg_V(t))))
}

real matrix mm_areg_V(struct mm_areg_struct scalar t)
{
    if (t.V==.z) t.V = mm_areg_XXinv(t)*mm_areg_s(t)^2
    return(t.V)
}

real matrix mm_areg_XXinv(struct mm_areg_struct scalar t)
{
    if (t.XXinv==.z) {
        t.XXinv = mm_ls_XXinv(t.ls) \ -(mm_areg_means(t) * mm_ls_XXinv(t.ls))
        t.XXinv = (t.XXinv, (t.XXinv[t.k+1,]' \
                   1/mm_areg_N(t) - mm_areg_means(t) * t.XXinv[t.k+1,]'))
    }
    return(t.XXinv)
}

real scalar mm_areg_rss(struct mm_areg_struct scalar t)
{
    if (t.rss==.z) {
        t.rss = quadcross(*t.w, mm_areg_e(t):^2)
    }
    return(t.rss)
}

real scalar mm_areg_ymean(struct mm_areg_struct scalar t)
{
    real scalar n
    
    if (t.ymean==.z) {
        n = rows(*t.y)
        if (n==0)               t.ymean = .
        else if (rows(*t.w)==1) t.ymean = quadsum(*t.y) / n
        else                    t.ymean = quadsum(*t.w :* *t.y) / mm_areg_N(t)
    }
    return(t.ymean)
}

real rowvector mm_areg_means(struct mm_areg_struct scalar t)
{
    if (t.means==.z) {
        if      (t.k==0)        t.means = J(1, 0, .)
        else if (rows(*t.X)==0) t.means = J(1, t.k, .)
        else                    t.means = quadcross(*t.w, *t.X) / mm_areg_N(t)
    }
    return(t.means)
}

real matrix mm_areg_yd(struct mm_areg_struct scalar t)
{
    return(t.yd)
}

real matrix mm_areg_Xd(struct mm_areg_struct scalar t)
{
    return(t.Xd)
}

real colvector mm_areg_omit(struct mm_areg_struct scalar t)
{
    return(mm_ls_omit(t.ls) \ 0)
}

real colvector mm_areg_k_omit(struct mm_areg_struct scalar t)
{
    return(mm_ls_k_omit(t.ls))
}

real scalar mm_areg_N(struct mm_areg_struct scalar t)
{
    if (t.N==.z) {
        t.N = rows(*t.w)==1 ? *t.w * rows(*t.y) : quadsum(*t.w)
    }
    return(t.N)
}

real colvector mm_areg_levels(struct mm_areg_struct scalar t)
{
    return(t.g->levels)
}

real colvector mm_areg_k_levels(struct mm_areg_struct scalar t)
{
    return(rows(t.g->idx))
}

real colvector mm_areg_n(struct mm_areg_struct scalar t)
{
    return(t.g->n)
}

end


*! {smcl}
*! {marker mm_aqreg}{bf:mm_aqreg.mata}{asis}
*! version 1.0.1  18sep2021  Ben Jann
version 9.2
mata:

// Main functions

real matrix mm_aqregfit(real colvector y, real colvector id, | real matrix X,
    real colvector w, real vector tau, real scalar sort, real scalar qd)
{
    if (args()<3) X = .
    if (args()<4) w = 1
    if (args()<5) tau = 0.5
    return(mm_aqreg_b(mm_aqreg(y, id, X, w, tau, sort, qd)))
}

struct mm_aqreg_struct {
    pointer(real colvector) scalar y, w
    pointer(real matrix) scalar    X
    real scalar                    k, N, ymean
    real rowvector                 means
    real rowvector                 tau, q
    real scalar                    beta0, gamma0
    real colvector                 alpha, delta
    struct mm_ls_struct scalar     L, S
    pointer(struct mm_areg_struct_g scalar) scalar g
}

struct mm_aqreg_struct scalar mm_aqreg(
    real colvector   y,
    real colvector   id,
    | real matrix    X,
      real colvector w,
      real vector    tau,
      real scalar    sort,
      real scalar    qd)
{
    if (args()<3) X = .
    if (args()<4) w = 1
    if (args()<5) tau = 0.5
    if (rows(id)!=rows(y)) _error(3200)
    return(_mm_aqreg(_mm_areg_g(id, sort), y, X, w, tau, qd))
}

struct mm_aqreg_struct scalar _mm_aqreg(
    struct mm_areg_struct_g scalar g,
    real colvector                 y,
    | real matrix                  X,
      real colvector               w,
      real vector                  tau,
      real scalar                  qd)
{
    real colvector ym, yd, r, rm, rd, e, u
    real matrix    Xm, Xd
    struct mm_aqreg_struct scalar t
    
    // setup
    if (args()<3) X = .
    if (args()<4) w = 1
    if (args()<5) tau = 0.5
    t.N = t.ymean = t.means = .z
    t.g = &g
    t.y = &y
    if (X==.) t.X = &(J(rows(y), 0, .))
    else      t.X = &X
    t.w = &w
    t.k = cols(*t.X)
    
    // quantiles
    if (rows(tau)>1) t.tau = tau'
    else             t.tau = tau
    if (length(tau)==0) _error(3200)
    if (any(tau:<=0 :| tau:>=1)) _error(3300, "tau must be between 0 and 1")
    
    // transform data
    ym = y; _mm_areg_gmean(g, ym, w, qd); yd = y - ym
    if (t.k) {; Xm = *t.X; _mm_areg_gmean(g, Xm, w, qd); Xd = *t.X - Xm; }
    else        Xd = Xm = *t.X
    
    // location fit (beta0: global constant; alpha: fixed effects)
    t.L = mm_ls(yd, Xd, w, 0, qd)
    t.L.y = t.L.X = t.L.w = NULL // release pointers to data (save memory)
    if (t.k) t.beta0 = mm_aqreg_ymean(t) - mm_ls_xb(t.L, mm_aqreg_means(t))
    else     t.beta0 = mm_aqreg_ymean(t)
    if (t.k) t.alpha = ym - mm_ls_xb(t.L, Xm)
    else     t.alpha = ym
    t.alpha = t.alpha :- t.beta0
    if (t.k) e = yd - mm_ls_xb(t.L, Xd)
    else     e = yd
    
    // scale fit (gamma0: global constant; delta: fixed effects)
    r = 2 * e :* ((e:>=0) :- mean(e:>=0, w))
    rm = r; _mm_areg_gmean(g, rm, w, qd); rd = r - rm
    t.S = mm_ls(rd, Xd, w, 0, qd)
    t.S.y = t.S.X = t.S.w = NULL // release pointers to data (save memory)
    if (t.k) t.gamma0 = mean(r, w) - mm_ls_xb(t.S, mm_aqreg_means(t))
    else     t.gamma0 = mean(r, w)
    if (t.k) t.delta = rm - mm_ls_xb(t.S, Xm)
    else     t.delta = rm
    if (t.k) u = e :/ (mm_ls_xb(t.S, *t.X) + t.delta)
    else     u = e :/ t.delta
    t.delta = t.delta :- t.gamma0
    
    // quantiles
    if (missing(u)) t.q = mm_quantile(select(u, u:<.),
                            rows(w)!=1 ? select(w, u:<.) : w, t.tau)
    else            t.q = mm_quantile(editmissing(u,0), w, t.tau)
    
    // return
    return(t)
}

// Results

real matrix mm_aqreg_xb(struct mm_aqreg_struct scalar t, | real matrix X)
{
    real matrix b
    
    b = mm_aqreg_b(t)
    if (args()==2) {
        if (cols(X)!=t.k) _error(3200)
        if (t.k) return(X * b[|1,1 \ t.k,.|] :+ b[t.k+1,])
        return(J(rows(X), 1, b))
    }
    if (t.k) return(*t.X * b[|1,1 \ t.k,.|] :+ b[t.k+1,])
    return(J(rows(*t.y), 1, b))
}

real matrix mm_aqreg_ue(struct mm_aqreg_struct scalar t)
{
    return(*t.y :- mm_aqreg_xb(t))
}

real matrix mm_aqreg_xbu(struct mm_aqreg_struct scalar t)
{
    return(mm_aqreg_xb(t) + mm_aqreg_u(t))
}

real matrix mm_aqreg_u(struct mm_aqreg_struct scalar t)
{
    return(t.alpha :+ t.delta :* J(rows(t.delta), 1, t.q))
}

real matrix mm_aqreg_e(struct mm_aqreg_struct scalar t)
{
    return(*t.y :- mm_aqreg_xbu(t))
}

real scalar mm_aqreg_ymean(struct mm_aqreg_struct scalar t)
{
    if (t.ymean==.z) {
        if (rows(*t.y)==0) t.ymean = .
        else               t.ymean = quadcross(*t.w, *t.y) / mm_aqreg_N(t)
    }
    return(t.ymean)
}

real rowvector mm_aqreg_means(struct mm_aqreg_struct scalar t)
{
    if (t.means==.z) {
        if      (t.k==0)        t.means = J(1, 0, .)
        else if (rows(*t.X)==0) t.means = J(1, t.k, .)
        else                    t.means = quadcross(*t.w, *t.X) / mm_aqreg_N(t)
    }
    return(t.means)
}

real matrix mm_aqreg_b(struct mm_aqreg_struct scalar t)
{
    return((mm_ls_b(t.L) :+ mm_ls_b(t.S) :* J(t.k, 1, t.q)) \ 
           (t.beta0 :+ t.gamma0 :* t.q))
}

real colvector mm_aqreg_omit(struct mm_aqreg_struct scalar t)
{
    return(mm_ls_omit(t.L) \ 0)
}

real colvector mm_aqreg_k_omit(struct mm_aqreg_struct scalar t)
{
    return(mm_ls_k_omit(t.L))
}

real scalar mm_aqreg_N(struct mm_aqreg_struct scalar t)
{
    if (t.N==.z) {
        t.N = rows(*t.w)==1 ? *t.w * rows(*t.y) : quadsum(*t.w)
    }
    return(t.N)
}

real colvector mm_aqreg_levels(struct mm_aqreg_struct scalar t)
{
    return(t.g->levels)
}

real colvector mm_aqreg_k_levels(struct mm_aqreg_struct scalar t)
{
    return(rows(t.g->idx))
}

real colvector mm_aqreg_n(struct mm_aqreg_struct scalar t)
{
    return(t.g->n)
}

real rowvector mm_aqreg_tau(struct mm_aqreg_struct scalar t)
{
    return(t.tau)
}

real rowvector mm_aqreg_q(struct mm_aqreg_struct scalar t)
{
    return(t.q)
}

real colvector mm_aqreg_beta(struct mm_aqreg_struct scalar t)
{
    return(mm_ls_b(t.L) \ t.beta0)
}

real colvector mm_aqreg_alpha(struct mm_aqreg_struct scalar t)
{
    return(t.alpha)
}

real colvector mm_aqreg_gamma(struct mm_aqreg_struct scalar t)
{
    return(mm_ls_b(t.S) \ t.gamma0)
}

real colvector mm_aqreg_delta(struct mm_aqreg_struct scalar t)
{
    return(t.delta)
}

end


*! {smcl}
*! {marker mm_rbinomial}{bf:mm_rbinomial.mata}{asis}
*! version 1.0.1, Ben Jann, 14apr2006
version 9.1
mata:

real matrix mm_rbinomial(real matrix n, real matrix p)
{
	real matrix x
	real scalar r, R, c, C
	transmorphic scalar rn, cn, rp, cp

	R = max((rows(n),rows(p)))
	C = max((cols(n),cols(p)))
	rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r))
	cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c))
	rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r))
	cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c))
	x = J(R,C,.)
	for (r=1;r<=R;r++) {
		for (c=1;c<=C;c++) {
			x[r,c] = _mm_rbinomial(n[*rn,*cn], p[*rp,*cp])
		}
	}
	return(x)
}

real scalar _mm_rbinomial(real scalar n, real scalar p)
{
	if (p<0|p>1) return(.)
	if (n<=0|(n-trunc(n))!=0) return(.)
//rejection technique
	if (n<50|p>.03) {
		return(colsum(uniform(n,1):<p))
	}
//geometric distribution technique
	real scalar sum, x
	sum = 0
	for (x=1;x<=n;x++) {
		sum = sum + ceil((ln(uniform(1,1))/ln(1-p))-1)
		if (sum > n-x) break
	}
	return(x-1)
}

end


*! {smcl}
*! {marker mm_cebinomial}{bf:mm_cebinomial.mata}{asis}
*! version 1.0.1, Ben Jann, 14apr2006
version 9.1
mata:

real matrix mm_cebinomial(real matrix n, real matrix k, real matrix p)
{
	real matrix x
	real scalar r, R, c, C
	transmorphic scalar rn, cn, rk, ck, rp, cp

	R = max((rows(n),rows(k),rows(p)))
	C = max((cols(n),cols(k),cols(p)))
	rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r))
	cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c))
	rk = (rows(k)==1 ? &1 : (rows(k)<R ? _error(3200) : &r))
	ck = (cols(k)==1 ? &1 : (cols(k)<C ? _error(3200) : &c))
	rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r))
	cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c))
	x = J(R,C,.)
	for (r=1;r<=R;r++) {
		for (c=1;c<=C;c++) {
			x[r,c] = _mm_cebinomial(n[*rn,*cn], k[*rk,*ck], p[*rp,*cp])
		}
	}
	return(x)
}

real scalar _mm_cebinomial(real scalar n, real scalar k, real scalar p)
{
	real scalar e, e0, i

	if (p<0|p>1) return(.)                 //success prob. out of range
	if (n<=0|(n-trunc(n))!=0) return(.)    //n out of range or non-integer
	if (k<0|k>n|(k-trunc(k))!=0) return(.) //k out of range or non-integer
	if (k==0) return(n*p)
	e = 0
	e0 = 0
	for (i=1;i<=n-k;i++) {
		e = e + Binomial(n,k+i,p)
		if (e==e0) break
		e0 = e
	}
	return(e / Binomial(n,k,p) + k)
}
end


*! {smcl}
*! {marker mm_benford}{bf:mm_benford.mata}{asis}
*! version 1.0.0, Ben Jann, 19jun2007
version 9.2
mata:

real vector mm_benford(
 real vector digit0,
 | real scalar pos0,
   real scalar base0)
{
    real scalar    pos, base
    real vector    digit

// check arguments and set defaults
    if (args()<3 | base0>=.) base = 10
    else {
        base = trunc(base0)
        if (base <0 | base >10) _error(3300, "base must be in [2,10]")
    }
    if (args()<2 | pos0>=.) pos = 1
    else {
        pos = trunc(pos0)
        if (pos<1) _error(3300, "pos must be 1 or larger")
    }
    digit = trunc(digit0)
    if (any(digit:<0) | any(digit:>=base))
      _error(3300, "digit must be in [0,base-1] (base is "+strofreal(base)+")")

// case 1: pos==1 (simple formula)
    if (pos==1) return(ln(1 :+ 1:/digit)/ln(base))

// case 2: pos > 1
    if (cols(digit)==1) return( rowsum(ln(1 :+ 1:/(
     J(rows(digit),1,base)*(base^(pos-2)..base^(pos-1)-1) :+ digit))) / ln(base))
    else                return( colsum(ln(1 :+ 1:/(
     (base^(pos-2)::base^(pos-1)-1)*J(1,cols(digit),base) :+ digit))) / ln(base))
}

end


*! {smcl}
*! {marker mm_cauchy}{bf:mm_cauchy.mata}{asis}
version 9.2
// CFBaum 3aug2007
// minor changes by Ben Jann, 09aug2007
// from Wikipedia
// http://en.wikipedia.org/wiki/Cauchy_distribution
mata:
// density function of the Cauchy-Lorentz distribution
real matrix mm_cauchyden(real matrix x, real scalar x0, real scalar gamma)
{
    real matrix cauchyden

    if (gamma <= 0) {
        return(x*.)
    }
    cauchyden = 1 :/ ( pi() * gamma :* ( 1 :+ ( (x :- x0) :/ gamma ) :^2) )
    return(cauchyden)
}

// cumulative distribution function of the Cauchy-Lorentz distribution
real matrix mm_cauchy(real matrix x, real scalar x0, real scalar gamma)
{
    real matrix cauchy

    if (gamma <= 0) {
        return(x*.)
    }
    cauchy = 1 / pi() :* atan( (x :- x0) :/ gamma ) :+ 0.5
    return(cauchy)
}

// tail probability of the Cauchy-Lorentz distribution
// note: cauchy(x,0,1) is the standard Cauchy distribution
//       with cauchytail(x,0,1) = ttail(1,x)
real matrix mm_cauchytail(real matrix x, real scalar x0, real scalar gamma)
{
    return(1 :- mm_cauchy(x, x0, gamma))
}

// inverse cumulative distribution function of the Cauchy-Lorentz distribution
real matrix mm_invcauchy(real matrix p, real scalar x0, real scalar gamma)
{
    real matrix invcauchy

    if ( gamma <= 0 | min(p) < 0 | max(p) > 1 ) {
        return(p*.)
    }
    invcauchy = x0 :+ gamma :* tan( pi() :* ( p :- 0.5 ) ) :+ ln(p:>=0:&p:<=1)
    return(invcauchy)
}
end


*! {smcl}
*! {marker mm_subset}{bf:mm_subset.mata}{asis}
*! version 1.0.1, Ben Jann, 02may2007
version 9.2
mata:

/*-SUBSETS---------------------------------------------------------------*/
// Generate all combinations (subsets) one-by-one (lexicographic)
// based on Algorithm 5.8 from Reingold et al. (1977)
// Usage:
//
//        info = mm_subsetsetup(n,k)
//        while ((c = mm_subset(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// The total number of combinations can be computed using
// Mata's -comb()- function, i.e.
//
//        Ctot = comb(n, k)

struct mm_subsetinfo {
	real scalar    n, k, i, j, counter, algorithm
	real colvector x, y
}

struct mm_subsetinfo scalar mm_subsetsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s
// setup
	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
	else if (s.k>s.n) _error(3300,"k may not be larger than n")
	if (s.n<1|s.k<1) s.i = 0  // => return empty set
	else {
		s.x = -1 \ (1::s.k)
		s.i = 1
	}
	s.i = s.i + 1 // offset i by 1
	s.k = s.k + 1 // offset k by 1
	s.counter = 0
	return(s)
}

real colvector mm_subset(struct mm_subsetinfo scalar s)
{
	real scalar    j
	real colvector subset

	if (s.i==1) return(J(0,1,.)) /*done*/
// get current subset (to be returned at end)
	subset = s.x[|2 \ s.k|]
// generate next subset
	s.i = s.k
	while (s.x[s.i] == s.n-s.k+s.i) {
		s.i = s.i - 1
	}
	s.x[s.i] = s.x[s.i] + 1
	for (j=s.i+1;j<=s.k;j++) s.x[j] = s.x[j-1] + 1
/* the following would improve speed but does not work for some reason
	if (s.i<s.k) s.x[|s.i+1 \ s.k|] = s.x[|s.i \ s.k-1|] :+ 1
*/
	s.counter = s.counter + 1
	return(subset)
}

// wrapper for mm_subsetinfo()/mm_subset()
real matrix mm_subsets(real scalar n, | real scalar k)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	s = mm_subsetsetup(n, (k<. ? k : n))
	if ((set = mm_subset(s)) == J(0,1,.)) return(set)
	res = J(rows(set), comb(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_subset(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}


/*-COMPOSITIONS----------------------------------------------------------*/
// Generate all k-way compositions of n, one-by-one
// Default algorithm: direct (less elegant but faster) (anti-lexicographic)
// Alternate algorithm: indirect via subsets (lexicographic)
//
// Usage:
//
//        info = mm_compositionsetup(n,k)
//        while ((c = mm_composition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// or:
//
//        info = mm_compositionsetup(n,k,1)
//        while ((c = mm_composition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// mm_ncompositions() returns the total number of compositions

real scalar mm_ncompositions(
 real scalar n,
 | real scalar k0)
{
	real scalar k

	k = trunc(k0<. ? k0 :  n)
	return(comb(trunc(n)+k-1,k-1))
}

struct mm_subsetinfo scalar mm_compositionsetup(
 real scalar n,
 | real scalar k,
   real scalar alt) // !=0 => alternate algorithm
{
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	if (alt) {
		s = _mm_composition2setup(n,k)
		s.algorithm = 2
		return(s)
	}
	s = _mm_compositionsetup(n,k)
	s.algorithm = 1
	return(s)
}

real colvector mm_composition(struct mm_subsetinfo scalar s)
{
	if (s.algorithm==1) return(_mm_composition(s))
	if (s.algorithm==2) return(_mm_composition2(s))
	_error(3498, "unknown algorithm")
}

// default algorithm

struct mm_subsetinfo scalar _mm_compositionsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = (k<. ? trunc(k) :  s.n)
	if (s.n>=.) _error(3351)
	if (s.k<1) {
		s.i = 0
		s.x = J(0,1,.)
	}
	else if (s.n<1) {
		s.i = 0
		s.x = J(s.k,1,0)
	}
	else if (s.k<2) {
		s.i = 0
		s.x = J(s.k,1,s.n)
	}
	else {
		s.x = s.n \ J(s.k-1,1,0)
		s.i = 1
	}
	s.j = 2
	s.counter = 0
	return(s)
}

real colvector _mm_composition(struct mm_subsetinfo scalar s)
{
	real colvector subset

	subset = s.x
	if (s.i<1) {
		if (subset!=J(0,1,.)) s.counter = s.counter + 1
		s.x = J(0,1,.)
		return(subset) /*done*/
	}
	while (s.x[s.i]==0 & s.i>1) { // move back
		s.x[s.i] = s.x[s.j]
		s.x[s.j] = 0
		s.j = s.i--
	}
	if (s.x[s.i]>0) {
		s.x[s.i] = s.x[s.i]-1
		s.x[s.j] = s.x[s.j]+1
		if (s.j<s.k) s.i = s.j++
	}
	else {
		s.i = 0
		s.x = J(0,1,.)
	}
	s.counter = s.counter + 1
	return(subset)
}

// alternate algorithm

struct mm_subsetinfo scalar _mm_composition2setup(
 real scalar n,
 | real scalar k0)
{
	real scalar k

	k = trunc(k0<. ? k0 :  n)
	return(mm_subsetsetup(trunc(n)+k-1,k-1))
}

real colvector _mm_composition2(struct mm_subsetinfo scalar s)
{
	real colvector subset

	if (s.i==1) {
		if (s.k==1) {  // length of composition is 1
			s.k = 0
			s.counter = s.counter + 1
			return(s.n)
		}
		return(J(0,1,.)) /*done*/
	}
	subset = mm_subset(s)
	return( (subset \ s.n+1) - (0 \ subset) :- 1 )
}

// wrapper for mm_compositioninfo()/mm_composition()
real matrix mm_compositions(real scalar n, | real scalar k, real scalar alt)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	s = mm_compositionsetup(n, (k<. ? k : n), alt)
	if ((set = mm_composition(s)) == J(0,1,.)) return(set)
	res = J(rows(set), mm_ncompositions(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_composition(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}

/*-RANDOM SUBSETS--------------------------------------------------------*/
// Return random combination (subsets)
// Usage:
//
//        for (i=1;i<=reps;i++) {
//             c = mm_rsubset(n, k)
//             ... c ...
//        }
//
// mm_rsubset() returns jumbled sets; apply sort() if you need
// ordered sets, e.g.
//
//        c = sort(mm_rsubset(n, k), 1)

real colvector mm_rsubset(
 real scalar n,
 | real scalar k)
{
	if (k<. & k>n) _error(3300, "k may not be larger than n")
	if (n<1|k<1) return(J(0,1,.))
	return(mm_unorder2(n)[|1 \ (k<. ? k : n)|])
}

/*-RANDOM COMPOSITIONS---------------------------------------------------*/
// Return random combination (subsets) and random composition
// Usage:
//
//        for (i=1;i<=reps;i++) {
//             c = mm_rcomposition(n, k)
//             ... c ...
//        }

real colvector mm_rcomposition(
 real scalar n0,
 | real scalar k0)
{
	real scalar    n, k
	real colvector subset

	n = trunc(n0)
	k = (k0<. ? trunc(k0) : n)
	if (n>=.) _error(3351)
	if (k<1) return(J(0,1,.))
	if (n<1) return(J(k,1,0))
	if (k<2) return(J(k,1,n))
	subset = sort(mm_rsubset(n+k-1,k-1), 1)
	return( (subset \ n+k) - (0 \ subset) :- 1 )
}

/*-PARTITIONS------------------------------------------------------------*/
// Generate all k-way partitions of n, one-by-one
// Default algorithm: based on ZS1 algorithm in Zoghbi&Stojmenovic (1998)
//                    (anti-lexicographic)
// Alternate algorithm: based on Algorithm 5.12 in Reingold et al. (1977)
//                    (dictionary order)
// Usage:
//
//        info = mm_partitionsetup(n,k)
//        while ((c = mm_partition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// or:
//
//        info = mm_partitionsetup(n,k,1)
//        while ((c = mm_partition(info)) != J(0,1,.)) {
//             ... c ...
//        }
//
// mm_npartitions() returns the total number of compositions, based
// on algorithms from
// http://home.att.net/~numericana/answer/numbers.htm#partitions

// number of partitions of n with k or fewer addends (slow)
real scalar mm_npartitions(
 real scalar n,
 | real scalar k)
{
	real scalar      u, i, j, m
	real colvector   p, a

	m = min((n,k))
	if (m==n) return(_mm_npartitions(n))
	p = J(n+1,1,1)
	for (u=2;u<=m;u++) {
		a = p
		p = J(n+1,1,0)
		for (i=0;i<=n;i=i+u) {
			for (j=i+1;j<=n+1;j++) {
				p[j] = p[j] + a[j-i]
			}
		}
	}
	return(p[n+1])
}
// total number of partitions of n (fast)
real scalar _mm_npartitions(
 real scalar n)
{
	real scalar      i, j, s, k
	real colvector   p

	p = J(n+1,1,1)
	for (i=1;i<=n;i++) {
		j = 1 ; k = 1 ; s = 0
		while (j>0) {
			j = i - (3*k*k + k)/2
			if (j>=0) s = s - (-1)^k * p[j+1]
			j = i - (3*k*k - k)/2
			if (j>=0) s = s - (-1)^k * p[j+1]
			k = k + 1
		}
		p[i+1] = s
	}
	return(p[n+1])
}

struct mm_subsetinfo scalar mm_partitionsetup(
 real scalar n,
 | real scalar k,
   real scalar pad,  // zero-padded output
   real scalar alt)  // !=0 => alternate algorithm
{
	struct mm_subsetinfo scalar s

	if (args()<3) pad = 0
	if (args()<4) alt = 0
	if (alt) {
		s = _mm_partition2setup(n,k)
		s.algorithm = 3 + (pad!=0)
		return(s)
	}
	s = _mm_partitionsetup(n,k)
	s.algorithm = 1 + (pad!=0)
	return(s)
}

real colvector mm_partition(struct mm_subsetinfo scalar s)
{
	if (s.algorithm==1) return(_mm_partition(s,0))
	if (s.algorithm==2) return(_mm_partition(s,1))
	if (s.algorithm==3) return(_mm_partition2(s,0))
	if (s.algorithm==4) return(_mm_partition2(s,1))
	_error(3498, "unknown algorithm")
}

// default algorithm

struct mm_subsetinfo scalar _mm_partitionsetup(
 real scalar n,
 | real scalar k)
{
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
//	else if (s.k>s.n) _error(3300,"k may not be larger than n")
	if (s.n<1 | s.k<1 ) s.i = -1
	else {
		s.x = s.n \ J(s.k-1, 1, 1)
		s.i = 1
	}
	s.j = 1
	s.counter = 0
	return(s)
}

real colvector _mm_partition(        // original ZS1 Algorithm does not
 struct mm_subsetinfo scalar s,      // support k-way partitions; changes
 real scalar pad)                    // are indicated
{
	real scalar    r, t, l
	real colvector subset

	if (s.i==-1) return(J(0,1,.)) /*done*/
	subset = s.x[|1 \ s.j|]
	if (pad) subset = subset \ J(s.k-rows(subset),1,0)
	l = (s.k<s.n ? s.k : s.n)     // added
	if (s.j>=l) {                 // added; original stopping rule is (s.x[1]==1)
		if (s.x[1]<=s.x[l]+1) {      // changed
			s.i = -1
			s.counter = s.counter + 1
			return(subset)
		}
		while (s.x[s.i]<=s.x[l]+1) {  // step back loop added
			s.j = s.j + s.x[s.i] - 1
			s.i = s.i - 1
		}
	}
	if (s.x[s.i]==2) {
		s.j      = s.j + 1
		s.x[s.i] = 1
		s.x[s.j] = 1             // added
		s.i      = s.i - 1
	}
	else {
		r = s.x[s.i] - 1
		t = s.j - s.i + 1
		s.x[s.i] = r
		while (t>=r) {
			s.i = s.i + 1
			s.x[s.i] = r
			t = t - r
		}
		if (t==0) {
			s.j = s.i
		}
		else {
			s.j = s.i + 1
			if (t==1) s.x[s.i+1] = 1    // added
			if (t>1) {
				s.i    = s.i + 1
				s.x[s.i] = t
			}
		}
	}
	s.counter = s.counter + 1
	return(subset)
}

// alternate algorithm

struct mm_subsetinfo scalar _mm_partition2setup(
 real scalar n,
 | real scalar k)
{
	real scalar offset
	struct mm_subsetinfo scalar s

	s.n = trunc(n)
	s.k = trunc(k)
	if (s.n>=.) _error(3351)
	if (s.k>=.) s.k = s.n
	offset = 2
	if (s.n<1 | s.k<1 ) s.i = 0 + offset
	else {
		s.i = 1 + offset
		s.x = s.y = J(s.n+offset, 1, .)
		s.x[-1+offset] = 0
		s.y[-1+offset] = 0
		s.x[0+offset]  = 0
		s.y[0+offset]  = s.n + 1
		s.x[1+offset]  = s.n
		s.y[1+offset]  = 1
	}
	s.counter = 0
	return(s)
}

real colvector _mm_partition2(
 struct mm_subsetinfo scalar s,
 real scalar pad)
{
	real scalar    skip, sum
	real colvector subset

	skip = 1
	while (skip) {
		if (s.i<=2) return(J(0,1,.)) /*done*/
		if ((skip = (sum(s.x[|3 \ s.i|])>s.k))==0)
		 subset = _mm_partition2_expand(s.x[|3 \ s.i|], s.y[|3 \ s.i|],
		 (pad ? s.k : .))
		sum = s.x[s.i]*s.y[s.i]
		if (s.x[s.i]==1) {
			s.i = s.i - 1
			sum = sum + s.x[s.i]*s.y[s.i]
		}
		if (s.y[s.i-1]==s.y[s.i]+1) {
			s.i      = s.i - 1
			s.x[s.i] = s.x[s.i] + 1
		}
		else {
			s.y[s.i] = s.y[s.i] + 1
			s.x[s.i] = 1
		}
		if (sum>s.y[s.i]) {
			s.y[s.i+1] = 1
			s.x[s.i+1] = sum - s.y[s.i]
			s.i        = s.i + 1
		}
	}
	s.counter = s.counter + 1
	return(subset)
}

real colvector _mm_partition2_expand(
 real vector m,
 real vector p,
 real scalar k)
{
	real scalar    i, j, l
	real colvector res

	if (k<.) res = J(k,1,0)
	else     res = J(sum(m),1,0)
	l = 0
	for (i=1;i<=rows(p);i++) {
		for (j=1;j<=m[i];j++) {
			res[++l] = p[i]
		}
	}
	return(res)
}

// wrapper for mm_partitioninfo()/mm_partition()
real matrix mm_partitions(real scalar n, | real scalar k, real scalar alt)
{
	real scalar    i
	real colvector set
	real matrix    res
	struct mm_subsetinfo scalar s

	if (args()<3) alt = 0
	s = mm_partitionsetup(n, (k<. ? k : n), 1, alt)
	if ((set = mm_partition(s)) == J(0,1,.)) return(set)
	res = J(rows(set), mm_npartitions(n, rows(set)), .)
	res[,i=1] = set
	while ((set = mm_partition(s)) != J(0,1,.)) res[,++i] = set
	return(res)
}

end


*! {smcl}
*! {marker mm_mgof}{bf:mm_mgof.mata}{asis}
*! version 1.0.7  Ben Jann  12jan2008
*  mm_mgof: goodness-of-fit tests for multinomial data

version 9.2
mata:

real matrix mm_mgof(
    real colvector f,    // observed count
  | real colvector h0,   // expected distribution
    string scalar m0,    // method: "approx" (default), "mc", "rc", or "ee"
    string vector s0,    // test-stats: "x2", "lr", "cr", "mlnp", "ksmirnov"
    real scalar lambda,  // lambda for Cressie-Read statistic
    real scalar arg1,    // nfit for "approx"; reps for "mc"
    real scalar dots,    // display dots with mc/ee
    real scalar arg2)    // will be replaced; used by mgof_ee to return counter
{
    real scalar     i, n, hsum
    string scalar   m, s
    string vector   stats
    real vector     p
    pointer(real colvector) scalar       h
    pointer(real scalar function) vector s1

// parsing and defaults
    if (args()<7) dots = 0
    m = mm_strexpand(m0, ("approx", "mc", "rc", "ee"), "approx", 0)
    stats = ("x2", "lr", "cr")
    if (m!="approx") stats = (stats, "mlnp", "ksmirnov")
    if (length(s0)<1) s = "x2"
    else {
        s = J(length(s0),1,"")
        for (i=1;i<=length(s0);i++) {
            s[i] = mm_strexpand(s0[i], stats, "x2", 1)
        }
    }
    if (m!="approx" & lambda<0) _error(3498, "lambda<0 not allowed with "+m)

// check data / prepare null distribution
    if (missing(f)) _error(3351, "f has missing values")
    if (any(f:<0)) _error(3498, "f has negative values")
    if (any(f:<0)) _error(3498, "f has negative values")
    if (m=="ee" | m=="rc" | (m=="mc" & anyof(s, "mlnp"))) {
        if (f!=trunc(f)) _error(3498, "non-integer f not allowed in this context")
    }
    if (anyof(stats,"cr") & lambda<0 & anyof(f,0))
     _error(3498, "some f are 0: lambda<0 not allowed")
    if (args()<2) h0 = 1
    if (rows(f)!=rows(h0) & rows(h0)!=1) _error(3200)
    if (missing(h0)) _error(3351, "h has missing values")
    if (any(h0:<=0)) _error(3498, "some h are negative or 0")
    n = colsum(f)
    hsum = colsum(h0)
    if (n!=hsum) h = &(h0*n/hsum)
    else h = &h0

// return (0,1) if n. of obs is 0 or n. of categories < 2
    if ( n==0 | rows(f)<2 ) return(J(rows(s),1,0),J(rows(s),1,1))

// gather statistics functions
    s1 = J(rows(s),1,NULL)
    for (i=1;i<=length(s);i++) {
        if      (s[i]=="x2")        s1[i] = &_mm_mgof_x2()
        else if (s[i]=="lr")        s1[i] = &_mm_mgof_lr()
        else if (s[i]=="cr")        s1[i] = &_mm_mgof_cr()
        else if (s[i]=="mlnp")      s1[i] = &_mm_mgof_mlnp()
        else if (s[i]=="ksmirnov")  s1[i] = &_mm_mgof_ksmirnov()
        else _error(3498,s[i]+ " invalid")
    }

// approximate chi2-test
    if (m=="approx") return(_mm_mgof_approx(s1, f, *h, lambda, arg1))

// Monte Carlo test
    if (m=="mc") return(_mm_mgof_mc(s1, f, *h, lambda, arg1, dots))

// exhaustive enumeration / random composition test
    if (s=="mlnp") {
        if (m=="rc") return(_mm_mgof_rc(NULL, f, *h))
        return(_mm_mgof_ee(NULL, f, *h))
    }
    p = mm_which(s:=="mlnp")
    if (rows(p)>0) {
        p = p[1]
        p = p \ (p>1 ? (1::p-1) : J(0,1,.)) \ (p<rows(s) ? (p+1::rows(s)) : J(0,1,.))
        s1 = s1[p[| 2 \ .|]]
        p = invorder(p)
    }
    else p = (2::length(s)+1)
    if (m=="rc") return(_mm_mgof_rc(s1, f, *h, lambda, arg1)[p,])
    return(_mm_mgof_ee(s1, f, *h, lambda, anyof(s,"ksmirnov"), dots, arg2)[p,])
}

// large sample chi2-approximation test for multinomial distributions
real matrix _mm_mgof_approx(
 pointer(real scalar function) vector s, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar nfit)   // n. of fitted parameters (default 0)
{
    real scalar                     n, i, k, nstat
    real colvector                  stat
    pointer(real colvector) scalar  h

// compute hypothetical distribution if missing
    n = colsum(f)
    k = rows(f)
    if (rows(h0)<=1) h = &(n :* J(k,1,1/k))
    else             h = &h0
// compute statistics and p-values
    nstat = length(s)
    stat = J(nstat,1,.)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h,lambda)
    return(stat, chi2tail(rows(f) - (nfit<. ? nfit : 0) - 1,stat))
}

// Monte Carlo exact test for discrete distributions
// (based on sampling from the hypothetical distribution)
real matrix _mm_mgof_mc(
 pointer(real scalar function) vector s, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar reps0,  // replications (default=10000)
   real scalar dots)   // display dots
{
    real scalar                     n, j, i, k, nstat, reps, uni
    real scalar                     doti, ndots
    real colvector                  fj, H, stat, p, hp
    pointer(real colvector) scalar  h

// set defaults and compute hypothetical distribution
    if (args()<5 | reps0>=.) reps = 10000
    else                     reps = reps0
    if (args()<6) dots = 0
    n = round(colsum(f)) // so that, e.g., 10 are sampled if n=9.9999...
    k = rows(f)
    if (uni = (rows(h0)==1)) {
        uni = 1
        h  = &(n :* J(k,1,1/k))
        H  = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        uni = mm_isconstant(h0)
        H  = mm_colrunsum(h0)
        H  = H / H[rows(H)]
        hp = (h0) / colsum(h0)
        h  = &h0
    }
// compute statistics and p-values
    nstat = length(s)
    stat = p = J(nstat,1,0)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (reps<=0) return(stat, J(nstat,1,.))
    if (dots) {
        printf("\n{txt}Percent completed ({res}%g{txt} replications)\n",reps)
        display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
        ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
        doti = 1
    }
    for (j=1; j<=reps; j++) {
        if (uni) fj = _mm_mgof_mc_usmpl(n, k) // = mm_srswr(n,rows(f),1)
        else     fj = _mm_mgof_mc_smpl(n, H) // = mm_upswr(n, h, 1)  (slower)
        for (i=1; i<=nstat; i++)  p[i] = p[i] +
          (stat[i] <= (*s[i])(fj,*h, lambda, n, hp, H))
        if (dots) {
            if (j*50 >= doti*ndots*reps) {
                printf(ndots*".")
                displayflush()
                doti++
            }
        }
    }
    if (dots) display("")
    return(stat, p/reps)
}
// Monte Carlo subroutine: sample from uniform multinomial
real colvector _mm_mgof_mc_usmpl(
 real scalar n, // sample size
 real scalar k) // number of categories
{
    real scalar     i
    real colvector  u, res

    u = ceil(uniform(n,1)*k)
    res = J(k,1,0)
    if (n>1.5^(k-2)) { // faster for small k relative to n
        for (i=1; i<=k; i++) res[i] = colsum(u:==i)
    }
    else {
        for (i=1;i<=n;i++) res[u[i]] = res[u[i]] + 1
    }
    return(res)
}
// Monte Carlo subroutine: sample from non-uniform multinomial
real colvector _mm_mgof_mc_smpl(
 real scalar n,    // sample size
 real colvector F) // distribution function (cumulative)
{
    real scalar     i, j
    real colvector  u, res

    if (n>2.3^(rows(F)-3)) {  // faster for small k relative to n
        u = uniform(n,1)
        res = J(rows(F),1,0)
        res[1] = sum(u:<=F[1])
        for (j=2;j<=rows(F);j++) {
            res[j] = sum(u:>F[j-1] :& u:<=F[j])
        }
        return(res)
    }
    u = sort(uniform(n,1),1)
    j=1
    res = J(rows(F),1,0)
    for (i=1;i<=n;i++) {
        while (u[i]>F[j]) j++
        res[j] = res[j]+1
    }
    return(res)
}

// Random composition exact test for discrete distributions
// -> undocumented; do not use this
//    the method is highly unreliable if the number of potential
//    compositions is large
real matrix _mm_mgof_rc(
 pointer(real scalar function) vector s0, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar reps0)  // replications (default=10000)
{
    real scalar                           n, i, j, k, nstat, reps, lnscale
    real colvector                        fj, H, stat, statj, p, jj, se, hp
    pointer(real colvector) scalar        h
    pointer(real scalar function) vector  s

// set defaults and compute hypothetical distribution
    if (args()<5 | reps0>=.) reps = 10000
    else                     reps = reps0
    if (s0==NULL) s = &_mm_mgof_mlnp()
    else {
        if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0
        else             s = &_mm_mgof_mlnp() , s0
    }
    n = colsum(f)
    k = rows(f)
    if (rows(h0)==1) {
        h = &(n :* J(k,1,1/k))
        H = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        h = &h0
        H = mm_colrunsum(*h)
        H = H / H[rows(H)]
        hp = (*h) / colsum(*h)
    }
// compute statistics and p-values
    nstat = length(s)
    stat = statj = p = jj = se = J(nstat,1,0)
    for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (reps<=0) return(stat, J(nstat,2,.))
    for (j=1; j<=reps; j++) {
        fj = mm_rcomposition(n, k)
        for (i=1; i<=nstat; i++) {
            statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
            if (stat[i] <= statj[i]) {
                jj[i] = jj[i] + 1
                p[i]  = p[i] + (exp(-statj[1])-p[i])/jj[i]     // mean updating
                se[i] = se[i] + (exp(-statj[1])^2-se[i])/jj[i] // mean updating
            }
        }
    }
    p = p:*jj; se = se:*jj
    // lnscale = log. of inverse sampling probability
    lnscale = lnfactorial(n+k-1) - lnfactorial(k-1) - lnfactorial(n) - ln(reps)
    se = exp( lnscale :+ ln(sqrt(reps*(1/(reps-1) * se - (p/reps):^2))) )
//  se = mm_ncompositions(n, k) * sqrt((1/(reps-1) * se - (p/reps):^2)/reps)
    p  = exp( lnscale :+ ln(p) )
//  p  = mm_ncompositions(n, k) / reps * p
    return(stat, p, se)
}

// exhaustive enumeration exact test (performs the exact
// multinomial g.o.f. test plus, optionally, exact tests
// for specified additional statistics)
// WARNING: use this test for very small samples only since
// computation time is linear to the number of compositions
// =comb(n+k-1,k-1), where n is the sample size and k is
// the number of categories.
// IMPORTANT EXCEPTION: when the null distribution is uniform,
// compution time is determined by the number of partitions
// which is magnitudes smaller than the number of compositions
real matrix _mm_mgof_ee(
 pointer(real scalar function) vector s0, // statistics to be tested
 real colvector f,     // observed count
 real colvector h0,    // expected count
 | real scalar lambda, // lambda for Cressie-Read statistic
   real scalar force,  // do not use partitions (for ksmirnov)
   real scalar dots,   // display dots
   real scalar counter) // will be replaced
{
    real scalar                           n, i, k, nstat, uni, w
    real scalar                           doti, ndots, reps
    real colvector                        fj, j, H, stat, statj, p, hp
    pointer(real colvector) scalar        h
    pointer(real scalar function) vector  s
    struct mm_subsetinfo scalar           info

// set defaults and compute hypothetical distribution
    if (args()<6) dots = 0
    if (s0==NULL) s = &_mm_mgof_mlnp()
    else {
        if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0
        else             s = &_mm_mgof_mlnp() , s0
    }
    n = colsum(f)
    k = rows(f)
    uni = (force==1 ? 0 : 1)
    if (rows(h0)==1) {
        h = &(n :* J(k,1,1/k))
        H = rangen(1/k,1,k)
        hp = J(k,1,1/k)
    }
    else {
        h = &h0
        if (mm_isconstant(*h)==0) uni = 0 // reset to 0 if h not constant
        H = mm_colrunsum(*h)
        H = H / H[rows(H)]
        hp = (*h) / colsum(*h)
    }
// compute statistics and p-values
    nstat = length(s)
    stat = statj = p = j = J(nstat,1,0)
    for (i=1; i<=nstat; i++)  stat[i] = (*s[i])(f,*h, lambda, n, hp, H)
    if (uni) {         // uniform null distribution: work with partitions
        if (dots) {
            reps = mm_npartitions(n,k)
            printf("\n{txt}Percent completed ({res}%g{txt} partitions)\n",reps)
            display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
            ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                    (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
            doti = 1
        }
        info = _mm_partitionsetup(n, k)
        while ((fj = _mm_partition(info,1)) != J(0,1,.)) {
            w = exp(lnfactorial(info.k) - quadcolsum(lnfactorial(_mm_panels(fj))))
            for (i=1; i<=nstat; i++) {
                statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
                if (stat[i] <= statj[i]) {
                    j[i] = j[i] + w // statj[1] = -ln(p) => p = exp(-statj[1])
                    p[i] = p[i] + (exp(-statj[1])-p[i])*w/j[i]  // mean updating
                }
            }
            if (dots) {
                if (info.counter*50 >= doti*ndots*reps) {
                    printf(ndots*".")
                    displayflush()
                    doti++
                }
            }
        }
        if (dots) display("")
        counter = info.counter
        return(stat, p:*j)
    }
    if (dots) {
        reps = mm_ncompositions(n,k)
        printf("\n{txt}Percent completed ({res}%g{txt} compositions)\n",reps)
        display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hline 5} 100")
        ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 :
                (reps>=5 ? 10 : (reps>=2 ? 25 : 50)))))
        doti = 1
    }
    info = _mm_compositionsetup(n, k)
    while ((fj = _mm_composition(info)) != J(0,1,.)) {
        for (i=1; i<=nstat; i++) {
            statj[i] = (*s[i])(fj,*h, lambda, n, hp, H)
            if (stat[i] <= statj[i]) {
                j[i] = j[i] + 1  // statj[1] = -ln(p) => p = exp(-statj[1])
                p[i] = p[i] + (exp(-statj[1])-p[i])/j[i]  // mean updating
            }
        }
        if (dots) {
            if (info.counter*50 >= doti*ndots*reps) {
                printf(ndots*".")
                displayflush()
                doti++
            }
        }
    }
    if (dots) display("")
    counter = info.counter
    return(stat, p:*j)
}

// Pearson's chi2 statistic subroutine
real scalar _mm_mgof_x2(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    return(quadcolsum((f-h):^2:/h))
}

// log likelihood ratio statistic subroutine
real scalar _mm_mgof_lr(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    return(2*quadcolsum(f:*(ln(f)-ln(h))))
}

// Cressie-Read statistic
real scalar _mm_mgof_cr(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | real scalar lambda,       // lambda parameter; default is 2/3
      transmorphic matrix opt2, // not used
      transmorphic matrix opt3, // not used
      transmorphic matrix opt4) // not used
{
    real scalar L
    L = (lambda<. ? lambda : 2/3)
    if (abs(L)<1e-6)        return( 2 * quadcolsum(f:*(ln(f)-ln(h))) )
    else if (abs(L+1)<1e-6) return( 2 * quadcolsum(h:*(ln(h)-ln(f))) )
    else if (L>0) return( 2/(L*(L+1)) * quadcolsum(f:*((f:/h):^L :- 1)) )
                  return( 2/(L*(L+1)) * quadcolsum(f:*((h:/f):^-L :- 1)) )
}

// -ln(p) statistic subroutine (multinomial test)
real scalar _mm_mgof_mlnp(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      real scalar n,            // n. of obs [ =colsum(f) ]
      real colvector p,         // hypothical propotions [ = h/colsum(h) ]
      transmorphic matrix opt4) // not used
{
    if (args()>=5) return(-(lnfactorial(n) -
     quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(p))))
    return(-(lnfactorial(quadcolsum(f)) -
     quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(h/quadcolsum(h)))))
}

// kolmogorov-smirnov statistic subroutine
real scalar _mm_mgof_ksmirnov(
    real colvector f,           // observed count
    real colvector h,           // expected count
    | transmorphic matrix opt1, // not used
      real scalar n,            // n. of obs [ =colsum(f) ]
      transmorphic matrix opt3, // not used
      real colvector H)         // hypothetical cumulative [ = runsum(h)/sum(h) ]
{
    if (args()>=6) return(max(abs(H-mm_colrunsum(f)/n)))
    return(max(abs(mm_colrunsum(h)/colsum(h)-mm_colrunsum(f)/colsum(f))))
}


end


*! {smcl}
*! {marker mm_root}{bf:mm_root.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006
*! translation of zeroin.c from http://www.netlib.org/c/brent.shar

* Quoted from http://www.netlib.org/c/brent.shar:
* X * Algorithm
* X *    G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
* X *    computations. M., Mir, 1980, p.180 of the Russian edition
* X *
* X *    The function makes use of the bissection procedure combined with
* X *    the linear or quadric inverse interpolation.
* X *    At every step program operates on three abscissae - a, b, and c.
* X *    b - the last and the best approximation to the root
* X *    a - the last but one approximation
* X *    c - the last but one or even earlier approximation than a that
* X *        1) |f(b)| <= |f(c)|
* X *        2) f(b) and f(c) have opposite signs, i.e. b and c confine
* X *           the root
* X *    At every step Zeroin selects one of the two new approximations, the
* X *    former being obtained by the bissection procedure and the latter
* X *    resulting in the interpolation (if a,b, and c are all different
* X *    the quadric interpolation is utilized, otherwise the linear one).
* X *    If the latter (i.e. obtained by the interpolation) point is
* X *    reasonable (i.e. lies within the current interval [b,c] not being
* X *    too close to the boundaries) it is accepted. The bissection result
* X *    is used in the other case. Therefore, the range of uncertainty is
* X *    ensured to be reduced at least by the factor 1.6

* Main changes compared to the original:
* 1. mm_root() returns an error code and stored the solution in x.
* 2. The maximum # of iteration may be set.
* 3. The tolerance argument is optional.
* 4. Additional arguments may be passed on to f.
* 5. Special action is taken if f returns missing, convergence is not
*    reached, or f(ax) and f(bx) do not have opposite signs.

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_root(
 transmorphic x,      // bj: will be replaced by solution
 pointer(real scalar function) scalar f,
                      // Address of the function whose zero will be sought for
 real scalar ax,      // Root will be sought for within a range [ax,bx]
 real scalar bx,      //
 | real scalar tol,   // Acceptable tolerance for the root value (default 0)
   real scalar maxit, // bj: maximum # of iterations (default: 1000)
   `opts')            // bj: additional args to pass on to f
{
    transmorphic  fs            // setup for f
    real scalar   a, b, c       // Abscissae, descr. see above
    real scalar   fa, fb, fc    // f(a), f(b), f(c)
    real scalar   prev_step     // Distance from the last but one
    real scalar   tol_act       // Actual tolerance
    real scalar   p             // Interpolation step is calcu-
    real scalar   q             // lated in the form p/q; divi-
                                // sion operations is delayed
                                // until the last moment
    real scalar   new_step      // Step at this iteration
    real scalar   t1, cb, t2
    real scalar   itr

    if (args()<5) tol = 0       // bj: set tolerance
    if (args()<6) maxit = 1000  // bj: maximum # of iterations

    fs = mm_callf_setup(f, args()-6, `opts') // bj: prepare function call

    x = .                       // bj: initialize output

    a = ax;  b = bx;  fa = mm_callf(fs, a);  fb = mm_callf(fs, b)
    c = a;  fc = fa

    if ( fa==. ) return(0)      // bj: abort if fa missing

    if ( (fa > 0 & fb > 0) |    // bj: f(ax) and f(bx) do not
         (fa < 0 & fb < 0) ) {  //     have opposite signs
        if ( abs(fa) < abs(fb) ) {
            x = a; return(2)    // bj: fa closer to zero than fb
        }
        x = b; return(3)        // bj: fb closer to zero than fa
    }

    for (itr=1; itr<=maxit; itr++) {
        if ( fb==. ) return(0)            // bj: abort if fb missing

        prev_step = b-a

        if( abs(fc) < abs(fb) ) {         // Swap data for b to be the
            a = b;  b = c;  c = a;        // best approximation
            fa = fb;  fb = fc;  fc = fa
        }

        tol_act = 2*epsilon(b) + tol/2
        new_step = (c-b)/2

        if( abs(new_step) <= tol_act | fb == 0 ) {
             x = b                        // Acceptable approx. is found
             return(0)
        }

        // Decide if the interpolation can be tried
        if( abs(prev_step) >= tol_act     // If prev_step was large enough
             & abs(fa) > abs(fb) ) {      // and was in true direction,
                                          // Interpolation may be tried
            cb = c-b
            if( a==c ) {                  // If we have only two distinct
                t1 = fb/fa                // points linear interpolation
                p = cb*t1                 // can only be applied
                q = 1.0 - t1
            }
            else {                        // Quadric inverse interpolation
                q = fa/fc;  t1 = fb/fc;  t2 = fb/fa
                p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) )
                q = (q-1.0) * (t1-1.0) * (t2-1.0)
            }
            if( p>0 )                     // p was calculated with the op-
              q = -q                      // posite sign; make p positive
            else                          // and assign possible minus to
              p = -p                      // q
            if( p < (0.75*cb*q-abs(tol_act*q)/2) // If b+p/q falls in [b,c]
               & p < abs(prev_step*q/2) ) // and isn't too large
             new_step = p/q               // it is accepted
                                          // If p/q is too large then the
                                          // bisection procedure can
                                          // reduce [b,c] range to more
                                          // extent
        }

        if( abs(new_step) < tol_act ) {   // Adjust the step to be not less
            if( new_step > 0 )            // than tolerance
             new_step = tol_act
            else
             new_step = -tol_act
        }

        a = b;  fa = fb                   // Save the previous approx.
        b = b + new_step;  fb = mm_callf(fs, b) // Do step to a new approxim.
        if( (fb > 0 & fc > 0) | (fb < 0 & fc < 0) )  {
                      // Adjust c for it to have a sign opposite to that of b
            c = a;  fc = fa
        }
    }
    x = b
    return(1)                             // bj: convergence not reached
}

end


*! {smcl}
*! {marker mm_nrroot}{bf:mm_nrroot.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_nrroot(
 real scalar x,              // initial guess; will be replaced by solution
 pointer(function) scalar f, // the function; must return f() and f'()
 | real scalar tol,          // acceptable tolerance (default 0)
   real scalar maxit,        // maximum # of iterations (default: 1000)
   `opts')                   // additional args to pass on to f
{
    transmorphic  fs         // setup for f
    real scalar   itr
    real scalar   fx         // fx = ( f(x), f'(x) )
    real scalar   dx         // dx = fx[1] / fx[2]
    real scalar   tol_act    // actual tolerance

    if (args()<3) tol = 0      // default tolerance
    if (args()<4) maxit = 1000 // default maximum # of iterations

    fs = mm_callf_setup(f, args()-4, `opts') // prepare function call
    for (itr=1; itr<=maxit; itr++) {
        tol_act = 2e+4*epsilon(x) + tol/2
        fx = mm_callf(fs, x)
        dx = fx[1] / fx[2]
        x = x - dx
        if ( abs(dx) <= tol_act | missing(x) ) return(0)
    }
    return(1)                // convergence not reached
}

end


*! {smcl}
*! {marker mm_minim}{bf:mm_minim.mata}{asis}
*! version 1.0.0  Ben Jann  05aug2020
*! translation of function Brent_fmin() from file optimize.c included in the
*! source of R 4.0.2

/*
Quoted from optimize.c:
    an approximation  x  to the point where  f  attains a minimum  on
the interval  (ax,bx)  is determined.

INPUT..

ax    left endpoint of initial interval
bx    right endpoint of initial interval
f     function which evaluates  f(x, info)  for any  x
      in the interval  (ax,bx)
tol   desired length of the interval of uncertainty of the final
      result ( >= 0.)

OUTPUT..

fmin  abcissa approximating the point where  f  attains a minimum

    The method used is a combination of  golden  section  search  and
successive parabolic interpolation.  convergence is never much slower
than  that  for  a  Fibonacci search.  If  f  has a continuous second
derivative which is positive at the minimum (which is not  at  ax  or
bx),  then  convergence  is  superlinear, and usually of the order of
about  1.324....
    The function  f  is never evaluated at two points closer together
than  eps*abs(fmin)+(tol/3), where eps is  approximately  the  square
root  of  the  relative  machine  precision.   if   f   is a unimodal
function and the computed values of   f   are  always  unimodal  when
separated  by  at least  eps*abs(x)+(tol/3), then  fmin  approximates
the abcissa of the global minimum of  f  on the interval  ax,bx  with
an error less than  3*eps*abs(fmin)+tol.  if   f   is  not  unimodal,
then fmin may approximate a local, but perhaps non-global, minimum to
the same accuracy.
    This function subprogram is a slightly modified  version  of  the
Algol  60 procedure  localmin  given in Richard Brent, Algorithms for
Minimization without Derivatives, Prentice-Hall, Inc. (1973).
*/

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

real scalar mm_minim(
    pointer(real scalar function) scalar f, // address of function to be minimized
    real scalar ax,   // minimum will be sought for within [ax,bx]
    real scalar bx,
  | real scalar tol0, // acceptable tolerance (length of uncertainty interval)
    `opts')           // additional args to pass on to f
{
    transmorphic  fs            // setup for f
    real scalar   a, b, c, d, e, p, q, r, u, v, w, x
    real scalar   t2, fu, fv, fw, fx, xm, eps, tol1, tol3, tol
    
    if (tol0<=0) _error(3300)
    tol = (tol0>=. ? epsilon(1)^0.25 : tol0) // default tolerance
    fs = mm_callf_setup(f, args()-4, `opts') // prepare function call
    c = (3 - sqrt(5)) * .5                   // squared inverse of the golden ratio
    eps = epsilon(1)
    tol1 = eps + 1 // the smallest 1.000... > 1
    eps = sqrt(eps)
    a = ax; b = bx
    v = w = x = a + c * (b - a)
    d = e = 0
    fv = fw = fx = mm_callf(fs, x)
    tol3 = tol / 3
    while (1) {
        xm = (a + b) * .5
        tol1 = eps * abs(x) + tol3
        t2 = tol1 * 2
        /* check stopping criterion */
        if (abs(x - xm) <= t2 - (b - a) * .5) break
        p = q = r = 0
        /* fit parabola */
        if (abs(e) > tol1) {
            r = (x - w) * (fx - fv)
            q = (x - v) * (fx - fw)
            p = (x - v) * q - (x - w) * r
            q = (q - r) * 2
            if (q > 0) p = -p
            else       q = -q
            r = e
            e = d
        }
        /* a golden-section step */
        if (abs(p) >= abs(q * .5 * r) | p <= q * (a - x) | p >= q * (b - x)) {
            if (x < xm) e = b - x
            else        e = a - x
            d = c * e
        }
        /* a parabolic-interpolation step */
        else {
            d = p / q
            u = x + d
            /* f must not be evaluated too close to ax or bx */
            if (u - a < t2 | b - u < t2) {
                d = tol1
                if (x >= xm) d = -d
            }
        }
        /* f must not be evaluated too close to x */
        if (abs(d) >= tol1) u = x + d
        else if (d > 0)     u = x + tol1
        else                u = x - tol1
        fu = mm_callf(fs, u)
        /*  update  a, b, v, w, and x */
        if (fu <= fx) {
            if (u < x) b = x
            else       a = x
            v = w;    w = x;   x = u
            fv = fw; fw = fx; fx = fu
        }
        else {
            if (u < x) a = u
            else       b = u
            if (fu <= fw | w == x) {
                v = w; fv = fw
                w = u; fw = fu
            } 
            else if (fu <= fv | v == x | v == w) {
                v = u; fv = fu
            }
        }
    }
    return(x)
}

end


*! {smcl}
*! {marker mm_finvert}{bf:mm_finvert.mata}{asis}
*! version 1.0.1, Ben Jann, 16may2014

version 9.2

mata:

real vector mm_finvert(
    y,
    f,     //     nr:      brent:
    o1,    //     fd       lo
    | o2,  //     x0       up
      tol,
      maxit,
      opt)
{
    if (args()<5) tol = 0
    if (args()<6) maxit = 1000
    if (ispointer(o1)) {
        if (args()<4) o2 = y   // set initial guess to y
        if (args()<7)
            return(_mm_finvert_nr(y, f, o1, o2, tol, maxit))
        else
            return(_mm_finvert_nr(y, f, o1, o2, tol, maxit, opt))
    }
    if (args()==3)
        _error(3001, "expected 4 to 6 arguments but received 3")
    if (args()<7)
        return(_mm_finvert_brent(y, f, o1, o2, tol, maxit))
    else 
        return(_mm_finvert_brent(y, f, o1, o2, tol, maxit, opt))
}

real vector _mm_finvert_nr(
    real vector y,
    pointer(real scalar function) scalar f,
    pointer(real scalar function) scalar df,
    real vector x,
    real scalar tol,
    real scalar maxit,
    | opt)
{
    real scalar     i, I, resi, rc
    real vector     res
    pointer scalar  ix

    I = length(y)
    ix = (length(x)==1 ? &1 : (length(x)<I ? _error(3200) : &i))

    res = J(rows(y), cols(y), .)
    for (i=1; i<=I; i++) {
        if (args()<7)
            rc = mm_nrroot(resi=x[*ix], &_mm_finvert_nr_fdf(),
                tol, maxit, y[i], f, df)
        else
            rc = mm_nrroot(resi=x[*ix], &_mm_finvert_nr_fdf(),
                tol, maxit, y[i], f, df, opt)
        if (rc)
         _error(3360, "failure to converge for element " + strofreal(i))
        res[i] = resi
    }
    return(res)
}

real rowvector _mm_finvert_nr_fdf(
    real scalar x,
    real scalar y,
    pointer(real scalar function) scalar f,
    pointer(real scalar function) scalar df,
    | opt)
{
    if (args()<5)
        return( (*f)(x) - y, (*df)(x) )
    else
        return( (*f)(x, opt) - y, (*df)(x, opt) )
}

real vector _mm_finvert_brent(
    real vector y,
    pointer(real scalar function) scalar f,
    real vector lo,
    real vector up,
    real scalar tol,
    real scalar maxit,
    | opt)
{
    real scalar     i, I, resi, rc
    real vector     res
    pointer scalar  il, iu

    I = length(y)
    il = (length(lo)==1 ? &1 : (length(lo)<I ? _error(3200) : &i))
    iu = (length(up)==1 ? &1 : (length(up)<I ? _error(3200) : &i))

    res = J(rows(y), cols(y), .)
    for (i=1; i<=I; i++) {
        if (args()<7)
            rc = mm_root(resi=., &_mm_finvert_brent_f(),
                lo[*il], up[*iu], tol, maxit, y[i], f)
        else
            rc = mm_root(resi=., &_mm_finvert_brent_f(),
                lo[*il], up[*iu], tol, maxit, y[i], f, opt)
        if (rc==1)
         _error(3360, "failure to converge for element " + strofreal(i))
        if (rc)
         _error(3498, "invalid search range for element " + strofreal(i))
        res[i] = resi
    }
    return(res)

}

real scalar _mm_finvert_brent_f(
    real scalar x,
    real scalar y,
    pointer(real scalar function) scalar f,
    | opt)
{
    if (args()<4)
        return( (*f)(x) - y )
    else
        return( (*f)(x, opt) - y )
}

end


*! {smcl}
*! {marker mm_integrate}{bf:mm_integrate.mata}{asis}
*! version 1.0.0  Ben Jann  26jan2014

version 9.2

local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"

mata:

// Simpson's rule (quadratic interpolation)
real colvector mm_integrate_sr(
    pointer(real scalar function) scalar f, // function to be integrated
    real scalar a,      // lower limit
    real scalar b,      // upper limit
    real scalar n,      // integration intervals
    | real scalar el,   // element wise evaluation
    `opts')             // additional args to pass on to f
{
    transmorphic    fs            // setup for f
    real scalar     d, i
    real colvector  x, w
    
    // check options
    if (missing(n))       _error(3300, "n: missing not allowed")
    if ((n<2) | mod(n,2)) _error(3300, "n: must be a positive multiple of two")
    if (missing(a))       _error(3300, "a: missing not allowed")
    if (missing(b))       _error(3300, "b: missing not allowed")
    if (a>=b)             _error(3300, "a must be strictly smaller than b")
    
    // prepare function call
    fs = mm_callf_setup(f, args()-5, `opts')
    
    // integration setup
    d = (b-a)/n
    x = rangen(a, b, n+1)
    w = 1 \ colshape((J(n/2,1,4), J(n/2,1,2)), 1)
    w[n+1] = 1
    
    // pass column vector
    if (el==0 | args()<5) {
        return(d/3 * quadcolsum(mm_callf(fs, x) :* w))
    }
    // element wise evaluation
    for (i=1; i<=(n+1); i++) {
        x[i] = mm_callf(fs, x[i])
    }
    return(d/3 * quadcolsum(x :* w))
}

// Simpson's 3/8 rule (cubic interpolation)
real colvector mm_integrate_38(
    pointer(real scalar function) scalar f, // function to be integrated
    real scalar a,      // lower limit
    real scalar b,      // upper limit
    real scalar n,      // integration intervals
    | real scalar el,   // element wise evaluation
    `opts')             // additional args to pass on to f
{
    transmorphic    fs            // setup for f
    real scalar     d, i
    real colvector  x, w
    
    // check options
    if (missing(n))       _error(3300, "n: missing not allowed")
    if ((n<3) | mod(n,3)) _error(3300, "n: must be a positive multiple of three")
    if (missing(a))       _error(3300, "a: missing not allowed")
    if (missing(b))       _error(3300, "b: missing not allowed")
    if (a>=b)             _error(3300, "a must be strictly smaller than b")
    
    // prepare function call
    fs = mm_callf_setup(f, args()-5, `opts')
    
    // integration setup
    d = (b-a)/n
    x = rangen(a, b, n+1)
    w = 1 \ colshape((J(n/3,1,3), J(n/3,1,3), J(n/3,1,2)), 1)
    w[n+1] = 1
    
    // pass column vector
    if (el==0 | args()<5) {
        return(3*d/8 * quadcolsum(mm_callf(fs, x) :* w))
    }
    // element wise evaluation
    for (i=1; i<=(n+1); i++) {
        x[i] = mm_callf(fs, x[i])
    }
    return(3*d/8 * quadcolsum(x :* w))
}

end


*! {smcl}
*! {marker mm_ipolate}{bf:mm_ipolate.mata}{asis}
*! version 2.0.0  08jul2020  Ben Jann
version 9.2
mata:

real colvector mm_ipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{
    real colvector p0, p1, y1
    
    if (args()<4) outer = 0
    if (rows(y0)!=rows(x0)) _error(3200)
    p0 = order(x0,1)
    p1 = order(x1,1)
    y1 = J(rows(x1),1,.)
    y1[p1] = _mm_ipolate(x0[p0], y0[p0], x1[p1], outer)
    return(y1)
}

real colvector _mm_ipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{   // x0 and x1 assumed sorted
    real scalar    n, i, ai, bi
    real colvector a, b, x, y

    if (args()<4) outer = 0
    n = rows(x0)
    if (rows(y0)!=n) _error(3200)
    if (n<=1) return(mm_fastipolate(x0, y0, x1, outer))
    // collapse data
    a = x0:!=(x0[n]\x0[|1\n-1|]) // tag first obs within panel
    a[1] = 1                     // (should x0 be constant)
    a = select(1::n, a)          // get index of first obs within panel
    if ((i=rows(a))==n) return(mm_fastipolate(x0, y0, x1, outer)) // no ties
    b = x0:!=(x0[|2\.|]\x0[1])   // tag last obs within panel
    b[n] = 1                     // (should x0 be constant)
    b = select(1::n, b)          // get index of last obs within panel
    x = y = J(i,1,.)
    for (;i;i--) {
        ai = a[i]; bi = b[i]
        x[i] = x0[ai]
        if (ai==bi) {
            y[i] = y0[ai]
            continue
        }
        y[i] = mean(y0[|ai\bi|],1)
    }
    return(mm_fastipolate(x, y, x1, outer))
}

real colvector mm_fastipolate(real colvector x0, real colvector y0,
    real colvector x1, | real scalar outer)
{   // x0 and x1 assumed sorted; x0 assumed unique
    real scalar    n, i, j
    real colvector y1

    if (args()<4) outer = 0
    j = rows(x0); n = i = rows(x1)
    if (rows(y0)!=j) _error(3200)
    if (j<1) return(J(n, 1, .))
    y1 = J(n,1,.)
    // handle values of x1 > max(x0)
    for (;i;i--) {
        if (x1[i]<=x0[j]) break
    }
    if (outer) {
        if (i<n) y1[|i+1\.|] = J(n-i, 1, y0[j])
    }
    // handle values of x1 within minmax(x0)
    for (;i;i--) {
        for (;j;j--) {
            if (x0[j]==x1[i]) {
                y1[i] = y0[j]
                break
            }
            if (x0[j]<x1[i]) {
                y1[i] = y0[j] + (y0[j+1]-y0[j]) * (x1[i]-x0[j])/(x0[j+1]-x0[j])
                break
            }
        }
        if (!j) break
    }
    // handle values of x1 < min(x0)
    if (outer) {
        if (i) y1[|1\i|] = J(i, 1, y0[1])
    }
    return(y1)
}

end


*! {smcl}
*! {marker mm_polint}{bf:mm_polint.mata}{asis}
*! version 1.0.0, Ben Jann, 12jul2006
version 9.2
mata:

real vector mm_polint(
 real vector xa,
 real vector ya,
 real vector x,
 | real scalar degree) // degree
{
	real scalar  n, m, i, j, a, b
	real vector  y

	if ( args()<4 ) degree = 1
	n = length(xa)
	if (length(ya)!=n) _error(3200)
	if ( degree<1 ) _error(3498, "degree must be 1 or larger")
	if ( degree>=n ) _error(3498, "degree must be smaller than length of data")
	if ( degree!=trunc(degree) ) _error(3498, "degree must be integer")

	m = degree + 1
	y = J(rows(x), cols(x), .)
	j = 0
	for (i=1; i<=length(x); i++) {
		mm_hunt(xa, x[i], j)
		a = min((max((trunc(j-(m-1)/2+.5), 1)), n+1-m))
		b = a+m-1
		y[i] = _mm_polint(xa[|a \ b|], ya[|a \ b|], x[i])
	}
	return(y)
}

real scalar _mm_polint(
// translation of -polint- from Press et al. (1992), Numerical
// Recipes in C, p. 109-110, http://www.numerical-recipes.com
 real vector xa,
 real vector ya,
 real scalar x)
{
	real scalar  i, m, n, ns
	real scalar  y, den, dif, dift, ho, hp, w
	real vector  c, d

	n = length(xa)
	if ( length(ya) != n ) _error(3200)

	ns = 1
	dif = abs(x-xa[1])
	c = J(1, n, .)
	d = J(1, n, .)
	for (i=1; i<=n; i++) {
		if ( (dift=abs(x-xa[i])) < dif) {
			ns  = i
			dif = dift
		}
		c[i] = ya[i]
		d[i] = ya[i]
	}
	y = ya[ns--]
	for (m=1; m<n; m++) {
		for (i=1; i<=n-m; i++) {
			ho = xa[i] - x
			hp = xa[i+m] - x
			w = c[i+1] - d[i]
			if ( (den=ho-hp) == 0.0)
			 _error(3498, "invalid input data") // This error can occur
			 // only if two input xa's are (to within roundoff) identical.
			den = w / den
			d[i] = hp * den
			c[i] = ho * den
		}
		if ( 2*ns < (n-m) ) y = y + c[ns+1]
		else y = y + d[ns--]
	}
	return(y)
}

end


*! {smcl}
*! {marker mm_sqrt}{bf:mm_sqrt.mata}{asis}
*! version 1.0.0  29may2017  Christopher Baum and Ben Jann
version 9.2
mata:

real matrix mm_sqrt(real matrix A)
{
    real colvector  L
    real matrix     X
    pragma unset    L
    pragma unset    X
    
    if (isfleeting(A)) _symeigensystem(A, X, L)
    else                symeigensystem(A, X, L)
    return(X * diag(sqrt(L)) * X')
}

end


*! {smcl}
*! {marker mm_plot}{bf:mm_plot.mata}{asis}
*! version 1.0.0, Ben Jann, 03aug2006
version 9.2
mata:

void mm_plot(
 real matrix X,
 | string scalar type,
   string scalar opts)
{
	real scalar rc

	rc = _mm_plot(X, type, opts)
	if (rc) _error(rc)
}

real scalar _mm_plot(
 real matrix X,
 | string scalar type,
   string scalar opts)
{
	real scalar   n, N, k, j, rc, updata
	string scalar vars, plottype
	real vector   vid

	updata = st_updata()
	plottype = type
	if (plottype=="") plottype = "scatter"
	n = rows(X)
	k = cols(X)
	N = st_nobs()
	if (N<n & k>0) st_addobs(n-N)
	vid = J(1, k, .)
	for (j=1; j<=k; j++) {
		st_store((1,n), vid[j]=st_addvar("double", st_tempname()), X[,j])
		st_varlabel(vid[j], "Column " + strofreal(j))
		vars = vars + " " + st_varname(vid[j])
	}
	st_updata(updata)
	rc = _stata("twoway " + plottype + vars + ", " + opts)
	updata = st_updata() // twoway might change data-have-changed flag
	if (N<n & k>0) st_dropobsin((N+1,n))
	st_dropvar(vid)
	st_updata(updata)
	return(rc)
}

end


*! {smcl}
*! {marker mm_group}{bf:mm_group.mata}{asis}
*! version 1.0.0  09jul2020  Ben Jann
version 9.2
mata:

real colvector mm_group(transmorphic matrix X, | real rowvector idx)
{
    real scalar    r, c
    real colvector p
    
    r = rows(X); c = cols(X)
    if (args()<2) {
        if (r<=1 | c==0) return(J(r,1,1))
        p = order(X, 1..c)
        p[p] = _mm_group(X[p,])
    }
    else {
        if (length(idx)>c) _error(3200)
        if (r<=1 | c==0) return(J(r,1,1))
        p = order(X, idx)
        p[p] = _mm_group(X[p, abs(idx)])
    }
    return(p)
}

real colvector _mm_group(transmorphic matrix X)
{
    real scalar r, c
    real colvector p
    
    r = rows(X); c = cols(X)
    if (r<=1 | c==0) return(J(r,1,1))
    if (c==1) p = (X :!= (X[r] \ X[|1\r-1|]))
    else      p = (rowsum(X :!= (X[r,] \ X[|1,1 \ r-1,.|])) :!= 0)
    return(mm_colrunsum(p))
}

end




*! {smcl}
*! {marker mm_panels}{bf:mm_panels.mata}{asis}
*! version 1.0.3, Ben Jann, 03may2007
version 9.2
mata:

void function mm_panels(
 transmorphic vector X1,       // upper level ID variable (e.g. strata)
 transmorphic matrix info1,    // will be replaced (unless X1 and X2 missing)
 | transmorphic vector X2,     // lower level ID variable (e.g. clusters)
   transmorphic matrix info2)  // will be replaced (unless X2 missing)
{
	real scalar i, j, nc, b, e, b2, e2

	if (length(X1)>0 & X1!=.) info1 = _mm_panels(X1)
	if (args()<3) return
	if (length(X2)==0 | X2==.) {
		if (length(X1)>0 & X1!=.) info1 = info1, info1
		return
	}
	if (length(X1)==0 | X1==.) {
		info2 = _mm_panels(X2)
		info1 = colsum(info2), rows(info2)
		return
	}
	if (length(X1)!=length(X2)) _error(3200)
	info1 = info1, J(rows(info1),1,.)
	e = 0
	for (i=1; i<=rows(info1); i++) {
		nc = 1
		b = e + 2
		e = e + info1[i,1]
		for (j=b; j<=e; j++) {
			if (X2[j]!=X2[j-1]) nc++
		}
		info1[i,2] = nc
	}
	if (args()<4) return
	info2 = J(colsum(info1[.,2]), 1, .)
	e = e2 = 0
	for (i=1; i<=rows(info1); i++) {
		b = e + 1
		e = e + info1[i,1]
		b2 = e2 + 1
		e2 = e2 + info1[i,2]
		info2[|b2 \ e2|] = _mm_panels(X2[|b \ e|], info1[i,2])
	}
}

real colvector _mm_panels(transmorphic vector X, | real scalar np)
{
	real scalar i, j, r
	real colvector res

	r = length(X)
	if (r<1) return(J(0,1,.))
	if (args()<2) np = r
	res = J(np, 1, 1)
	j = 1
	for (i=2; i<=r; i++) {
		if (X[i]==X[i-1]) res[j] = res[j] + 1
		else j++
	}
	if (j==r) return(res)
	return(res[|1 \ j|])
}

// /*Old (slow) version of _mm_panels()*/
//real colvector _mm_panels(transmorphic vector X, | real scalar np)
//{
//	real scalar i, j, n
//	real colvector res
//
//	if (args()<2) np = mm_npanels(X)
//	if (length(X)==0) return(J(0,1,.))
//	res = J(np, 1, .)
//	n = j = 1
//	for (i=2; i<=length(X); i++) {
//		if (X[i]!=X[i-1]) {
//			res[j++] = n
//			n = 1
//		}
//		else n++
//	}
//	res[j] = n
//	return(res)
//}

real scalar mm_npanels(vector X)
{
	real scalar i, np

	if (length(X)==0) return(0)
	np = 1
	for (i=2; i<=length(X); i++) {
		if (X[i]!=X[i-1]) np++
	}
	return(np)
}

end


*! {smcl}
*! {marker mm_nunique}{bf:mm_nunique.mata}{asis}
*! version 1.1.1  09oct2020  Ben Jann
version 9.2
mata:

real scalar mm_nunique(transmorphic vector X)
{
    return(sum(mm_unique_tag(X)))
}

transmorphic vector mm_unique(transmorphic vector X, | real scalar order)
{
    if (rows(X)*cols(X)==0)  return(X)
    if (order==1 | order==2) return(select(X, mm_unique_tag(X, order)))
    // sorted
    if (cols(X)==1)    return(sort(select(X, mm_unique_tag(X)),1))
    if (!iscomplex(X)) return(sort(select(X, mm_unique_tag(X))',1)')
    return(transposeonly(sort(transposeonly(select(X, mm_unique_tag(X))),1)))
}

/*
real vector mm_unique_idx(transmorphic vector X, | real scalar order)
{
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (cols(X)==1) return(select(1::rows(X), mm_unique_tag(X, order)))
                    return(select(1..cols(X), mm_unique_tag(X, order)))
}
*/

real vector mm_unique_tag(transmorphic vector X, | real scalar order)
{
    real vector p
    
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (order==1) {
        // use stable sort order and tag first occurrence
        if (cols(X)==1)        p = mm_order(X,1,1)
        else if (iscomplex(X)) p = mm_order(transposeonly(X),1,1)'
        else                   p = mm_order(X',1,1)'
        p[p] = _mm_unique_tag(X[p])
    }
    else if (order==2) {
        // use stable sort order and tag last occurrence
        if (cols(X)==1)        p = mm_order(X,1,1)
        else if (iscomplex(X)) p = mm_order(transposeonly(X),1,1)'
        else                   p = mm_order(X',1,1)'
        p[p] = _mm_unique_tag(X[p], 1)
    }
    else {
        // use non-stable sort order (tag random element within ties)
        if (cols(X)==1)        p = order(X,1)
        else if (iscomplex(X)) p = order(transposeonly(X),1)'
        else                   p = order(X',1)'
        p[p] = _mm_unique_tag(X[p])
    }
    return(p)
}

real scalar _mm_nunique(transmorphic vector X)
{
    return(sum(_mm_unique_tag(X)))
}

transmorphic vector _mm_unique(transmorphic vector X)
{
    if (rows(X)*cols(X)==0) return(X)
    return(select(X, _mm_unique_tag(X)))
}

/*
real vector _mm_unique_idx(transmorphic vector X, | real scalar last)
{
    if (args()<2) last = 0
    if (rows(X)*cols(X)==0) return(J(rows(X),cols(X),.))
    if (cols(X)==1) return(select(1::rows(X), _mm_unique_tag(X, last)))
                    return(select(1..cols(X), _mm_unique_tag(X, last)))
}
*/

real vector _mm_unique_tag(transmorphic vector X, | real scalar last)
{
    real scalar r, c
    real vector p

    if (args()<2) last = 0
    r = rows(X); c = cols(X)
    if (r*c==0)      return(J(r,c,.))
    if (r==1 & c==1) return(J(1,1,1))
    if (last) {
        if (c==1) {
            p = (X :!= (X[|2\.|] \ X[1]))
            p[r] = 1 // should X be constant
            return(p)
        }
        p = (X :!= (X[|2\.|] , X[1]))
        p[c] = 1 // should X be constant
        return(p)
    }
    if (c==1) {
        p = (X :!= (X[r] \ X[|1\r-1|]))
        p[1] = 1 // should X be constant
        return(p)
    }
    p = (X :!= (X[c] , X[|1\c-1|]))
    p[1] = 1 // should X be constant
    return(p)
}

real scalar mm_nuniqrows(transmorphic matrix X)
{
    return(sum(mm_uniqrows_tag(X)))
}

transmorphic matrix mm_uniqrows(transmorphic matrix X, | real scalar order)
{
    if (rows(X)==0) return(X)
    if (order==1 | order==2) return(select(X, mm_uniqrows_tag(X, order)))
    return(mm_sort(select(X, mm_uniqrows_tag(X))))
}

/*
real vector mm_uniqrows_idx(transmorphic matrix X, | real scalar order)
{
    if (rows(X)==0) return(J(0,1,.))
    return(select(1::rows(X), mm_uniqrows_tag(X, order)))
}
*/

real vector mm_uniqrows_tag(transmorphic matrix X, | real scalar order)
{
    real vector p
    
    if (rows(X)==0) return(J(0,1,.))
    if (order==1) {
        p = mm_order(X,.,1)
        p[p] = _mm_uniqrows_tag(X[p,])
    }
    else if (order==2) {
        p = mm_order(X,.,1)
        p[p] = _mm_uniqrows_tag(X[p,], 1)
    }
    else {
        p = mm_order(X,.,0)
        p[p] = _mm_uniqrows_tag(X[p,])
    }
    return(p)
}

real scalar _mm_nuniqrows(transmorphic matrix X)
{
    return(sum(_mm_uniqrows_tag(X)))
}

transmorphic matrix _mm_uniqrows(transmorphic matrix X)
{
    if (rows(X)==0) return(X)
    return(select(X, _mm_uniqrows_tag(X)))
}

/*
real vector _mm_uniqrows_idx(transmorphic matrix X, | real scalar last)
{
    if (args()<2) last = 0
    if (rows(X)==0) return(J(0,1,.))
    return(select(1::rows(X), _mm_uniqrows_tag(X, last)))
}
*/

real vector _mm_uniqrows_tag(transmorphic matrix X, | real scalar last)
{
    real scalar r, c
    real vector p

    if (args()<2) last = 0
    r = rows(X); c = cols(X)
    if (r==0) return(J(0,1,.))
    if (r==1) return(J(1,1,1))
    if (last) {
        if (c==0) p = J(r,1,0)
        else      p = (rowsum(X :!= (X[|2,1 \ .,.|] \ X[1,])):!=0)
        p[r] = 1 // should X be constant
        return(p)
    }
    if (c==0) p = J(r,1,0)
    else      p = (rowsum(X :!= (X[r,] \ X[|1,1 \ r-1,.|])):!=0)
    p[1] = 1 // should X be constant
    return(p)
}

end



*! {smcl}
*! {marker mm_diff}{bf:mm_diff.mata}{asis}
*! version 1.0.1  Ben Jann  11may2021
version 9.2
mata:

numeric vector mm_diff(numeric vector x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = length(x)
    if (n<=lag) {
        if (cols(x)!=1) return(J(1,0,missingof(x)))
        return(J(0,1,missingof(x)))
    }
    return(x[|1+lag \ .|] - x[|1 \ n-lag|])
}

numeric matrix mm_rowdiff(numeric matrix x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = cols(x)
    if (n<=lag)     return(J(rows(x),0,missingof(x)))
    if (rows(x)==0) return(J(0, n-lag, missingof(x)))
    return(x[|1,1+lag \ .,.|] - x[|1,1 \ .,n-lag|])
}

numeric matrix mm_coldiff(numeric matrix x, | real scalar lag0)
{
    real scalar n, lag
    
    if (lag0>=.) lag =  1
    else         lag = abs(trunc(lag0))
    n = rows(x)
    if (n<=lag)     return(J(0,cols(x),missingof(x)))
    if (cols(x)==0) return(J(n-lag, 0, missingof(x)))
    return(x[|1+lag,1 \ .,.|] - x[|1,1 \ n-lag,.|])
}

end


*! {smcl}
*! {marker mm_isconstant}{bf:mm_isconstant.mata}{asis}
*! version 1.0.2, Ben Jann, 16jul2020
version 9.0
mata:

real scalar mm_isconstant(transmorphic matrix X)
{
    if (length(X)<2) return(1)
    return(allof(X, X[1,1]))
}

real scalar mm_issorted(transmorphic vector x, | real scalar descending)
{
    if (ispointer(x)) _error(3250)
    if (length(x)<=1) return(1)
    if (args()<2) descending = 0
    if (descending) {
        if (isstring(x)) {
            if (cols(x)==1) return(all(x[|1 \ length(x)-1|] :>= x[|2 \ .|]))
                            return(all(x[|1 \ length(x)-1|] :>= x[|2 \ .|]))
        }
        if (cols(x)==1) return(all((.z \ x[|1 \ length(x)-1|]) :>= x))
                        return(all((.z, x[|1 \ length(x)-1|]) :>= x))
    }
    if (isstring(x)) {
        if (cols(x)==1) return(all(x[|1 \ length(x)-1|] :<= x[|2 \ .|]))
                        return(all(x[|1 \ length(x)-1|] :<= x[|2 \ .|]))
    }
    if (cols(x)==1) return(all(x :<= (x[|2 \ .|] \ .z)))
                    return(all(x :<= (x[|2 \ .|], .z)))
}

end


*! {smcl}
*! {marker mm_colrunsum}{bf:mm_colrunsum.mata}{asis}
*! version 1.0.9  09jul2020  Ben Jann
version 9.0
mata:

numeric matrix mm_colrunsum(numeric matrix A, | real scalar mis, real scalar qd)
{
    numeric matrix B

    if (args()<2) mis = 0
    if (stataversion()>=1000) {
        if (args()<3) qd = 0
        return(_mm_colrunsum_goto10(A, mis, qd))
    }
    if (isfleeting(A)) {
        _mm_colrunsum(A, mis)
        return(A)
    }
    _mm_colrunsum(B=A, mis)
    return(B)
}

void _mm_colrunsum(numeric matrix Z, real scalar mis)
{
    real scalar i

    if (mis==0) _editmissing(Z, 0)
    for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + Z[i,]

//    // running sum using mean-update formula; see Gould 2006, SJ 6(4)
//    if (mis==0) _editmissing(Z, 0)
//    if (rows(Z)<2) return
//    for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + (Z[i,]-Z[i-1,])/i
//    Z = Z :* (1::rows(Z))

}

numeric matrix _mm_colrunsum_goto10(numeric matrix Z, real scalar mis, real scalar qd)
{
    if (qd) return(_mm_quadcolrunsum10(Z, mis))
    return(_mm_colrunsum10(Z, mis))
}


end


*! {smcl}
*! {marker mm_linbin}{bf:mm_linbin.mata}{asis}
*! version 1.0.9  12aug2020  Ben Jann
version 9.0
mata:

real colvector mm_linbin(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real colvector  p

    p = mm_order(x, ., 1) // stable sort order
    if (rows(w)==1) return(__mm_linbin(x[p], g)*w)
                    return(__mm_linbin_w(x[p], w[p], g))
}

real colvector _mm_linbin(
    real colvector x,  // data points (assumed sorted)
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    if (rows(w)==1) return(__mm_linbin(x, g)*w)
                    return(__mm_linbin_w(x, w, g))
}

real colvector __mm_linbin(
    real colvector x,  // data points (assumed sorted)
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     i, j, g1, g0, i1, i0
    real colvector  c, xx

    j = rows(g)
    i = rows(x)
    if (i<1) return(J(j, 1, 0))
    if (j<2) return(i)
    c = J(j, 1, 0)
    // data above grid
    g0 = g[j]
    i1 = i
    for (;i;i--) {
        if (x[i]<g0) break
    }
    i0 = i+1
    c[j] = (i1-i0+1)
    // data within grid range
    g1 = g0
    for (j--; j; j--) {
        g1 = g0
        g0 = g[j]
        i1 = i
        for (;i;i--) {
            if (x[i]<g0) break
        }
        i0 = i+1
        if (i1<i0) continue // no data in current bin
        xx = x[|i0\i1|]
        c[j] =            sum(g1 :- xx) / (g1 - g0)
        c[j+1] = c[j+1] + sum(xx :- g0) / (g1 - g0)
    }
    // data below grid
    if (i) c[1] = c[1] + i
    return(c)
}

real colvector __mm_linbin_w(
    real colvector x,  // data points (assumed sorted)
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     i, j, g1, g0, i1, i0
    real colvector  c, xx, ww

    j = rows(g)
    i = rows(x)
    if (rows(w)!=i) _error(3200)
    if (i<1) return(J(j, 1, 0))
    if (j<2) return(sum(w))
    c = J(j, 1, 0)
    // data above grid
    g0 = g[j]
    i1 = i
    for (;i;i--) {
        if (x[i]<g0) break
    }
    i0 = i+1
    if (i0<=i1) c[j] = sum(w[|i0 \ i1|])
    // data within grid range
    g1 = g0
    for (j--; j; j--) {
        g1 = g0
        g0 = g[j]
        i1 = i
        for (;i;i--) {
            if (x[i]<g0) break
        }
        i0 = i+1
        if (i1<i0) continue // no data in current bin
        xx = x[|i0\i1|]
        ww = w[|i0\i1|]
        c[j] =            sum(ww:*(g1 :- xx)) / (g1 - g0)
        c[j+1] = c[j+1] + sum(ww:*(xx :- g0)) / (g1 - g0)
    }
    // data below grid
    if (i) c[1] = c[1] + sum(w[|1 \ i|])
    return(c)
}

real colvector mm_fastlinbin(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     g1, i, j, d, z, ng
    real colvector  c

    ng = rows(g)
    i  = rows(x)
    if (ng<1) _error(3200)
    if (i<1)  return(J(ng, 1, 0))
    if (ng<2) return(mm_nobs(x,w))
    if (rows(w)!=1) return(_mm_fastlinbin_w(x, w, g))
    g1 = g[1]
    d = (g[ng] - g1)/(ng-1)
    c = J(ng, 1, 0)
    for (; i; i--) {
        z = (x[i] - g1) / d
        if (z<=0) {
            c[1] = c[1] +  1
            continue
        }
        if (z>=(ng-1)) {
            c[ng] = c[ng] +  1
            continue
        }
        j = trunc(z) + 1
        c[j] = c[j] + (j-z)
        j++
        c[j] = c[j] + (2-j+z)
    }
    return(c*w)
}

real colvector _mm_fastlinbin_w(
    real colvector x,  // data points
    real colvector w,  // weights
    real colvector g)  // grid points (assumed sorted)
{
    real scalar     g1, i, j, d, z, ng
    real colvector  c

    ng = rows(g)
    i  = rows(x)
    if (rows(w)!=i) _error(3200)
    g1 = g[1]
    d = (g[ng] - g1)/(ng-1)
    c = J(ng, 1, 0)
    for (; i; i--) {
        z = (x[i] - g1) / d
        if (z<=0) {
            c[1] = c[1] +  w[i]
            continue
        }
        if (z>=(ng-1)) {
            c[ng] = c[ng] +  w[i]
            continue
        }
        j = trunc(z) + 1
        c[j] = c[j] + w[i] * (j-z)
        j++
        c[j] = c[j] + w[i] * (2-j+z)
    }
    return(c)
}

end


*! {smcl}
*! {marker mm_exactbin}{bf:mm_exactbin.mata}{asis}
*! version 1.0.6  12aug2020  Ben Jann
version 9.0
mata:

real colvector mm_exactbin(
    real colvector x,   // data
    real colvector w,   // weights
    real colvector g,   // grid (assumed sorted)
    | real scalar dir,  // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    real colvector  p

    if (args()<4) dir = 0
    if (args()<5) include = 0
    p = mm_order(x, ., 1) // stable sort order
    if (rows(w)==1) return(__mm_exactbin(x[p], g, dir, include)*w)
                    return(__mm_exactbin_w(x[p], w[p], g, dir, include))
}

real colvector _mm_exactbin(
    real colvector x,   // data (assumed sorted)
    real colvector w,   // weights
    real colvector g,   // grid (assumed sorted)
    | real scalar dir,  // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    if (args()<4) dir = 0
    if (args()<5) include = 0
    if (rows(w)==1) return(__mm_exactbin(x, g, dir, include)*w)
                    return(__mm_exactbin_w(x, w, g, dir, include))
}

real colvector __mm_exactbin(
    real colvector x,    // data (assumed sorted)
    real colvector g,    // grid (assumed sorted)
    real scalar dir,     // 0 left inclusive (default),  else right inclusive
    real scalar include) // include data outside grid in first/last bin
{
    real scalar     i, i1, j, gj
    real colvector  c

    j = rows(g)
    if (j<2) _error(3200)
    i = rows(x)
    if (i<1) return(J(j-1, 1, 0))
    if (include==0) {
        if (x[1]<g[1] | x[i]>g[j]) _error("data out of grid range")
    }
    if (j<3) return(i)
    j--
    c = J(j, 1, 0)
    if (dir==0) {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = i1 - i
        }
    }
    else {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<=gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = i1 - i
        }
    }
    if (i) c[1] = c[1] + i // data below grid
    return(c)
}

real colvector __mm_exactbin_w(
    real colvector x,    // data (assumed sorted)
    real colvector w,    // weights
    real colvector g,    // grid (assumed sorted)
    real scalar dir,     // 0 left inclusive (default),  else right inclusive
    real scalar include) // include data outside grid in first/last bin
{
    real scalar     i, i1, j, gj
    real colvector  c
    
    j = rows(g)
    if (j<2) _error(3200)
    i = rows(x)
    if (i<1) return(J(j-1, 1, 0))
    if (include==0) {
        if (x[1]<g[1] | x[i]>g[j]) _error("data out of grid range")
    }
    if (j<3) return(sum(w))
    j--
    c = J(j, 1, 0)
    if (dir==0) {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = sum(w[|i+1 \ i1|])
        }
    }
    else {
        for (; j; j--) {
            gj = g[j]
            i1 = i
            for (;i;i--) {
                if (x[i]<=gj) break
            }
            if (i1<=i) continue // no data in current bin
            c[j] = sum(w[|i+1 \ i1|])
        }
    }
    if (i) c[1] = c[1] + sum(w[|1 \ i|]) // data below grid
    return(c)
}

real colvector mm_fastexactbin(
    real colvector x,      // data points
    real colvector w,      // weights
    real colvector g,      // grid (assumed sorted)
    | real scalar dir,     // 0 left inclusive (default),  else right inclusive
      real scalar include) // include data outside grid in first/last bin
{
    real scalar     g1, i, j, n, d, xi
    real colvector  c
    real rowvector  minmax
    pointer scalar  wi

    if (args()<4) dir = 0
    if (args()<5) include = 0
    n = rows(g) - 1
    i = rows(x)
    if (i<1) return(J(n, 1, 0))
    if (include==0) {
        minmax = minmax(x)
        if (minmax[1]<g[1] | minmax[2]>g[n+1]) _error("data out of grid range")
    }
    if (n<2) return(mm_nobs(x,w))
    wi = (rows(w)==1 ? &1 : (rows(w)!=i ? _error(3200) : &i))
    g1 = g[1]
    d = (g[n+1] - g1)/n
    c = J(n, 1, 0)
    if (dir==0) {
        for (; i; i--) {
            xi = x[i]
            j = floor((xi - g1) / d) + 1
            // data below grid
            if (j<1) {
                c[1] = c[1] +  w[*wi]
                continue
            }
            // data above grid
            if (j>n) {
                c[n] = c[n] +  w[*wi]
                continue
            }
            // roundoff error 1: j erroneously rounded up
            if (xi<g[j]) {
                // not sure whether this can happen at all...
                if (j>1) c[j-1] = c[j-1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // roundoff error 2: j erroneously rounded down
            if (xi>=g[j+1]) {
                if (j<n) c[j+1] = c[j+1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // no roundoff error
            c[j] = c[j] + w[*wi]
        }
    }
    else {
        for (; i; i--) {
            xi = x[i]
            j = ceil((xi - g1) / d)
            // data below grid
            if (j<1) {
                c[1] = c[1] +  w[*wi]
                continue
            }
            // data above grid
            if (j>n) {
                c[n] = c[n] +  w[*wi]
                continue
            }
            // roundoff error 1: j erroneously rounded up
            if (xi<=g[j]) {
                if (j>1) c[j-1] = c[j-1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // roundoff error 2: j erroneously rounded down
            if (xi>g[j+1]) {
                // not sure whether this can happen at all...
                if (j<n) c[j+1] = c[j+1] + w[*wi]
                else     c[j]   = c[j] + w[*wi]
                continue
            }
            // no roundoff error
            c[j] = c[j] + w[*wi]
        }
    }
    return(c)
}

end


*! {smcl}
*! {marker mm_makegrid}{bf:mm_makegrid.mata}{asis}
*! version 1.0.4, Ben Jann, 26jun2006
version 9.0
mata:

real colvector mm_makegrid(
 real colvector x, // data points
 | real scalar M,  // number of bins (default: 401)
   real scalar e,  // extend range by e on each side
   real scalar min,  // min
   real scalar max)  // max
{
	real scalar a, b

	a = (min<. ? min : min(x) - (e<. ? e : 0))
	b = (max<. ? max : max(x) + (e<. ? e : 0))
	return( rangen(a, b, (M<. ? M : 512)) )
}

end


*! {smcl}
*! {marker mm_seq}{bf:mm_seq.mata}{asis}
*! version 1.0.0  Ben Jann  04aug2020

version 9.2
mata:

real colvector mm_seq(real scalar a, real scalar b, real scalar d)
{
    if (missing((a,b,d))) return(J(0,1,.))
    if (a<=b)             return( a :+ (0::(b-a)/abs(d)) * abs(d) )
                          return( a :- (0::(a-b)/abs(d)) * abs(d) )
}

end


*! {smcl}
*! {marker mm_cut}{bf:mm_cut.mata}{asis}
*! version 1.0.0, Ben Jann, 16mar2006
version 9.1
mata:

real colvector mm_cut(real colvector x, real vector at, | real scalar sorted)
{
	real scalar i, j
	real colvector result, p

	result = J(rows(x),1,.)
	j = length(at)
	if (args()<3) sorted = 0
	if (sorted) {
		for (i=rows(x); i>0; i--) {
			if (x[i]>=.) continue
			for (; j>0; j--) {
				if (at[j]<=x[i]) break
			}
			if (j>0) result[i,1] = at[j]
		}
		return(result)
	}
	p = order(x,1)
	for (i=rows(x); i>0; i--) {
		if (x[p[i]]>=.) continue
		for (; j>0; j--) {
			if (at[j]<=x[p[i]]) break
		}
		if (j>0) result[p[i],1] = at[j]
	}
	return(result)
}
end


*! {smcl}
*! {marker mm_posof}{bf:mm_posof.mata}{asis}
*! version 1.0.0, Ben Jann, 24mar2006
version 9.0
mata:

real scalar mm_posof(transmorphic vector haystack, transmorphic scalar needle)
{
	real scalar i

	for (i=1; i<=length(haystack); i++) {
		if (haystack[i]==needle) return(i)
	}
	return(0)
}
end


*! {smcl}
*! {marker mm_which}{bf:mm_which.mata}{asis}
*! version 1.0.2, Ben Jann, 17apr2007
version 9.2
mata:

real matrix mm_which(real vector I)
{
	if (cols(I)!=1) return(select(1..cols(I), I))
	else return(select(1::rows(I), I))
}

end


*! {smcl}
*! {marker mm_locate}{bf:mm_locate.mata}{asis}
*! version 1.0.0, Ben Jann, 07jul2006

version 9.2

mata:

void mm_locate(
// translation of -locate- from Press et al. (1992), Numerical
// Recipes in C, p. 117, http://www.numerical-recipes.com
// Description: "Given an array xx[1..n], and given a value x,
// returns a value j such that x is between xx[j] and xx[j+1].
// xx must be monotonic, either increasing or decreasing. j=0
// or j=n is returned to indicate that x is out of range."
 numeric vector xx,
 numeric scalar x,
 real scalar j)
{
	real scalar  ju, jm, jl
	real scalar  ascnd
	real scalar  n

	n = length(xx)
	jl = 0
	ju = n + 1
	ascnd = (xx[n] >= xx[1])
	while (ju-jl > 1) {
		jm = trunc((ju+jl)/2)
		if (x >= xx[jm] == ascnd)
		  jl = jm
		else
		  ju = jm
	}
	if (x == xx[1])     j = 1
	else if(x == xx[n]) j = n - 1
	else                j = jl
}

void mm_hunt(
// translation of -hunt- from Press et al. (1992), Numerical
// Recipes in C, p. 118-9, http://www.numerical-recipes.com
// Description: "Given an array xx[1..n], and given a value x,
// returns a value jlo such that x is between xx[jlo] and
// xx[jlo+1]. xx[1..n] must be monotonic, either increasing or
// decreasing. jlo=0 or jlo=n is returned to indicate that x is
// out of range. jlo on input is taken as the initial guess for
// jlo on output.
 numeric vector xx,
 numeric scalar x,
 real scalar jlo)
{
	real scalar  jm, jhi, inc
	real scalar  ascnd
	real scalar  n

	n = length(xx)
	ascnd=(xx[n] >= xx[1])
	if (jlo <= 0 | jlo > n) {
		jlo = 0
		jhi = n + 1
	}
	else {
		inc = 1
		if (x >= xx[jlo] == ascnd) {
			if (jlo == n) return
			jhi = jlo + 1
			while (x >= xx[jhi] == ascnd) {
				jlo = jhi
				inc = inc + inc
				jhi = jlo + inc
				if (jhi > n) {
					jhi = n + 1
					break
				}
			}
		}
		else {
			if (jlo == 1) {
				jlo = 0
				return
			}
			jhi = jlo--
			while (x < xx[jlo] == ascnd) {
				jhi = jlo
				inc = inc + inc
				if (inc >= jhi) {
					jlo = 0
					break
				}
				else  jlo = jhi - inc
			}
		}
	}
	while (jhi-jlo != 1) {
		jm = trunc((jhi+jlo)/2)
		if (x >= xx[jm] == ascnd)
		 jlo = jm
		else
		 jhi = jm
	}
	if (x == xx[n]) jlo = n - 1
	if (x == xx[1]) jlo = 1
}

end


*! {smcl}
*! {marker mm_clip}{bf:mm_clip.mata}{asis}
*! version 1.0.0  11jul2020  Ben Jann
version 9.2
mata:

real matrix mm_clip(real matrix X, real scalar min, real scalar max, | real scalar miss)
{
    real matrix R
    
    if (args()<4) miss = 0
    if (isfleeting(X)) {
        _mm_clip(X, min, max, miss)
        return(X)
    }
    R = X
    _mm_clip(R, min, max, miss)
    return(R)
}

void _mm_clip(real matrix X, real scalar min, real scalar max, | real scalar miss)
{
    real scalar i, j, r, x
    
    if (args()<4) miss = 0
    if (min>max) _error(3300)
    r = rows(X)
    if (miss) {
        for (j=cols(X); j; j--) {
            for (i=r; i; i--) {
                x = X[i,j]
                if (x>=.) continue
                if (x<min)      X[i,j] = min
                else if (x>max) X[i,j] = max
            }
        }
        return
    }
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            x = X[i,j]
            if (x<min)      X[i,j] = min
            else if (x>max) X[i,j] = max
        }
    }
}

real matrix mm_clipmin(real matrix X, real scalar min)
{
    real matrix R
    
    if (isfleeting(X)) {
        _mm_clipmin(X, min)
        return(X)
    }
    R = X
    _mm_clipmin(R, min)
    return(R)
}

void _mm_clipmin(real matrix X, real scalar min)
{
    real scalar i, j, r
    
    r = rows(X)
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            if (X[i,j]<min) X[i,j] = min
        }
    }
}

real matrix mm_clipmax(real matrix X, real scalar max, | real scalar miss)
{
    real matrix R
    
    if (args()<3) miss = 0
    if (isfleeting(X)) {
        _mm_clipmax(X, max, miss)
        return(X)
    }
    R = X
    _mm_clipmax(R, max, miss)
    return(R)
}

void _mm_clipmax(real matrix X, real scalar max, | real scalar miss)
{
    real scalar i, j, r, x
    
    if (args()<3) miss = 0
    r = rows(X)
    if (miss) {
        for (j=cols(X); j; j--) {
            for (i=r; i; i--) {
                x = X[i,j]
                if (x>=.) continue
                if (x>max) X[i,j] = max
            }
        }
        return
    }
    for (j=cols(X); j; j--) {
        for (i=r; i; i--) {
            if (X[i,j]>max) X[i,j] = max
        }
    }
}

end


*! {smcl}
*! {marker mm_cond}{bf:mm_cond.mata}{asis}
*! version 1.0.0  29feb2008  Ben Jann
version 9.2
mata:

transmorphic matrix mm_cond(real matrix x, transmorphic matrix y, transmorphic matrix z)
{
        transmorphic matrix res
        real scalar r, R, c, C
        transmorphic scalar rx, cx, ry, cy, rz, cz

        if (eltype(y) != eltype(z)) _error(3250)
        R = max((rows(x),rows(y),rows(z)))
        C = max((cols(x),cols(y),cols(z)))
        rx = (rows(x)==1 ? &1 : (rows(x)<R ? _error(3200) : &r))
        cx = (cols(x)==1 ? &1 : (cols(x)<C ? _error(3200) : &c))
        ry = (rows(y)==1 ? &1 : (rows(y)<R ? _error(3200) : &r))
        cy = (cols(y)==1 ? &1 : (cols(y)<C ? _error(3200) : &c))
        rz = (rows(z)==1 ? &1 : (rows(z)<R ? _error(3200) : &r))
        cz = (cols(z)==1 ? &1 : (cols(z)<C ? _error(3200) : &c))
        res = J(R,C, missingof(y))
        for (r=1;r<=R;r++) {
                for (c=1;c<=C;c++) {
                        res[r,c] = (x[*rx,*cx] ? y[*ry,*cy] : z[*rz,*cz])
                }
        }
        return(res)
}

end


*! {smcl}
*! {marker mm_expand}{bf:mm_expand.mata}{asis}
*! version 1.0.0, Ben Jann, 07jun2006
version 9.2
mata:

transmorphic matrix mm_expand(transmorphic matrix X, real vector wr,
 | real vector wc, real scalar sort)
{
	transmorphic matrix Xnew

	if (args()<3) wc = 1
	if (args()<4) sort = 0
	if (isfleeting(X)) {
		_mm_expand(X, wr, wc, sort)
		return(X)
	}
	else {
		_mm_expand(Xnew=X, wr, wc, sort)
		return(Xnew)
	}
}

void _mm_expand(transmorphic matrix X, real vector wr,
 real vector wc, real scalar sort)
{
	real scalar i, b, e, n, add
	real colvector f
	pointer scalar fi

// expand rows
	if (length(wr)==1 & sort==0) _mm_repeat(X, wr, 1)
	else {
		f = trunc(wr:-1) :* (wr:>0)
		_editmissing(f,0)
		n = rows(X)
		if (length(f)==1) {
			fi = &1
			add = n*f
		}
		else {
			if (n!=length(f)) _error(3200)
			fi = &i
			add = sum(f)
		}
		if (add>0) {
			X = X \ J(add,cols(X),missingof(X))
			if (cols(X)>0) {
				if (sort) f = f :+ 1
				b = rows(X)+1
				for (i=n; i>=1; i--) {
					if (f[*fi]<1) continue
					e = b - 1
					b = b - f[*fi]
					X[|b,1 \ e,.|] = X[J(f[*fi],1,i), .]
				}
			}
		}
	}
// expand columns
	if (length(wc)==1 & sort==0) {
		_mm_repeat(X, 1, wc)
		return
	}
	f = trunc(wc:-1) :* (wc:>0)
	_editmissing(f,0)
	n = cols(X)
	if (length(f)==1) {
		fi = &1
		add = n*f
	}
	else {
		if (n!=length(f)) _error(3200)
		fi = &i
		add = sum(f)
	}
	if (add>0) {
		X = X , J(rows(X),add,missingof(X))
		if (rows(X)>0) {
			if (sort) f = f :+ 1
			b = cols(X)+1
			for (i=n; i>=1; i--) {
				if (f[*fi]<1) continue
				e = b - 1
				b = b - f[*fi]
				X[|1,b \ .,e|] = X[., J(1,f[*fi],i)]
			}
		}
	}
}

transmorphic matrix mm_repeat(transmorphic matrix X, real scalar wr,
 | real scalar wc)
{
	transmorphic matrix Xnew

	if (args()<3) wc = 1
	if (isfleeting(X)) {
		_mm_repeat(X, wr, wc)
		return(X)
	}
	else {
		_mm_repeat(Xnew=X, wr, wc)
		return(Xnew)
	}
}

void _mm_repeat(transmorphic matrix X, real scalar wr,
 real scalar wc)
{
	real scalar i, rr, rc, r, c

	r = rows(X)
	c = cols(X)
	rr = max((trunc(wr-1),0))
	if (rr>0) {
		X = X \ J(rr*r, c, missingof(X))
		if (r>0 & c>0)
		 for (i=1; i<=rr; i++) X[|i*r+1,1 \ i*r+r,.|] = X[|1,1 \ r,.|]
	}
	rc = max((trunc(wc-1),0))
	if (rc>0) {
		X = X , J(rows(X), rc*c, missingof(X))
		if (r>0 & c>0)
		 for (i=1; i<=rc; i++) X[|1,i*c+1 \ .,i*c+c|] = X[|1,1 \ .,c|]
	}
}

end


*! {smcl}
*! {marker mm_sort}{bf:mm_sort.mata}{asis}
*! version 1.0.0  09jul2020  Ben Jann
version 9.2
mata:

transmorphic matrix mm_sort(transmorphic matrix X,
    | real rowvector idx, real scalar stable)
{
    if (args()<2) return(X[mm_order(X), .])
    if (args()<3) return(X[mm_order(X, idx), .])
    return(X[mm_order(X, idx, stable), .])
}

real colvector mm_order(transmorphic matrix X,
    | real rowvector idx, real scalar stable)
{
    real scalar r, c, l
    
    if (args()<3) stable = 0
    // redirect to order()
    c = cols(X)
    if (!stable) {
        if (args()>=2 & idx!=.) return(order(X, idx))
        return(order(X, (c>0 ? 1..c : J(1,0,.))))
    }
    // check idx
    l = (idx==. ? 0 : length(idx))
    if (l>c) _error(3200)
    // one row or less: nothing to do
    r = rows(X)
    if (r<=1) return(J(r,1,1))
    // zero columns: return index
    if (c==0) return(1::r)
    // X is real: append index for stable order
    if (isreal(X)) {
        if (l==0) return(order((X,(1::r)), 1..c+1))
        if (anyof(idx,0)) _error(3300)
        return(order((X[,abs(idx)],(1::r)), ((1..l):*sign(idx), l+1)))
    }
    // X is not real; need to create a group id
    if (l==0) return(order((mm_group(X), (1::r)), (1,2)))
    return(order((mm_group(X, idx), (1::r)), (1,2)))
}

end



*! {smcl}
*! {marker mm_unorder2}{bf:mm_unorder2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

real colvector mm_unorder2(real scalar n) return(order(uniform(n,2),(1,2)))

end


*! {smcl}
*! {marker mm_jumble2}{bf:mm_jumble2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

transmorphic matrix mm_jumble2(transmorphic matrix x)
{
	return(x[mm_unorder2(rows(x)), .])
}

end


*! {smcl}
*! {marker mm__jumble2}{bf:mm__jumble2.mata}{asis}
*! version 1.0.0, Ben Jann, 29mar2006
version 9.0
mata:

function mm__jumble2(transmorphic matrix x)
{
	_collate(x,mm_unorder2(rows(x)))
}

end


*! {smcl}
*! {marker mm_pieces}{bf:mm_pieces.mata}{asis}
*! version 1.0.2, Ben Jann, 01jun2015
version 9.2
mata:

string rowvector mm_pieces(
    string scalar s,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar         i, j, k, n, l, b, nobr
    string scalar c
    string rowvector    res

    nobr = ( args()>2 ? nobreak : 0 )
    if (stataversion()>=1400 & !(args()>3 ? ascii : 0)) {
        return(_mm_pieces_goto14(s, w, nobr))
    }
    l = strlen(s)
    if (l<2 | w>=l) return(strtrim(s))
    res = J(1, _mm_npieces(s, w, nobr, (args()>3 ? ascii : 0)), "")
    j = k = n = 0
    b = 1
    for (i=1; i<=l; i++) {
        c = substr(s, i, 1)
        if (j<1) { // skip to first nonblank character
            if (c==" ") {
                b++
                continue
            }
        }
        j++
        if (i==l) res[++n] = strtrim(substr(s, b, .))
        else {
            if (c==" ") k = i
            if (j>=w) {
                if (k<1) {
                    if (nobr) continue
                    k = i
                }
                else {
                    if (substr(s, i+1, 1)==" ") k = i
                }
                res[++n] = strtrim(substr(s, b, k-b+1))
                j = i - k
                b = k + 1
                k = 0
            }
        }
    }
    return(res)
}

real matrix mm_npieces(
    string matrix S,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar  i, j
    real matrix  res

    res = J(rows(S),cols(S),1)
    if (w>=.) return(res)
    for (i=1; i<=rows(S); i++) {
        for (j=1; j<=cols(S); j++) {
            res[i,j] = _mm_npieces(S[i,j], w, (args()>2 ? nobreak : 0), 
                                              (args()>3 ? ascii : 0))
        }
    }
    return(res)
}

real scalar _mm_npieces(
    string scalar s,
    real scalar w,
    | real scalar nobreak,
      real scalar ascii)
{
    real scalar   i, j, k, n, l, nobr
    string scalar c

    nobr = ( args()>2 ? nobreak : 0 )
    if (stataversion()>=1400 & !(args()>3 ? ascii : 0)) {
        return(_mm_npieces_goto14(s, w, nobr))
    }
    l = strlen(s)
    if (l<2 | w>=l) return(1)
    j = k = n = 0
    for (i=1; i<=l; i++) {
        c = substr(s, i, 1)
        if (j<1) { // skip to first nonblank character
            if (c==" ") continue
        }
        j++
        if (i==l) n++
        else {
            if (c==" ") k = i
            if (j>=w) {
                if (k<1) {
                    if (nobr) continue
                    k = i
                }
                else {
                    if (substr(s, i+1, 1)==" ") k = i
                }
                n++
                j = i - k
                k = 0
            }
        }
    }
    if (n==0) n++
    return(n)
}

string rowvector _mm_pieces_goto14(string scalar s,
    real scalar w,
    real scalar nobr)
{
    return(_mm_pieces14(s, w, nobr))
}

real scalar _mm_npieces_goto14(string scalar s,
    real scalar w,
    real scalar nobr)
{
    return(_mm_npieces14(s, w, nobr))
}

end


*! {smcl}
*! {marker mm_regexr}{bf:mm_regexr.mata}{asis}
*! version 1.0.0  31jan2017  Ben Jann
version 9.2
mata:

string matrix mm_regexr(string matrix x, string matrix y, string matrix z)
{
    string matrix       res
    real scalar         r, R, c, C
    transmorphic scalar rx, cx, ry, cy, rz, cz
    pragma unset        r
    pragma unset        c

    R = max((rows(x),rows(y),rows(z)))
    C = max((cols(x),cols(y),cols(z)))
    rx = (rows(x)==1 ? &1 : (rows(x)<R ? _error(3200) : &r))
    cx = (cols(x)==1 ? &1 : (cols(x)<C ? _error(3200) : &c))
    ry = (rows(y)==1 ? &1 : (rows(y)<R ? _error(3200) : &r))
    cy = (cols(y)==1 ? &1 : (cols(y)<C ? _error(3200) : &c))
    rz = (rows(z)==1 ? &1 : (rows(z)<R ? _error(3200) : &r))
    cz = (cols(z)==1 ? &1 : (cols(z)<C ? _error(3200) : &c))
    res = J(R,C, "")
    for (r=1;r<=R;r++) {
        for (c=1;c<=C;c++) {
            res[r,c] = _mm_regexr(x[*rx,*cx], y[*ry,*cy], z[*rz,*cz])
        }
    }
    return(res)
}

string scalar _mm_regexr(string scalar s0, string scalar from, string scalar to)
{
    real scalar         i, j
    string scalar       s, t, BSLASH
    string rowvector    sub
    pragma unset        s
    
    if (regexm(s0, from)) {
        sub = regexs()
        sub = (sub, J(1, 10-cols(sub), ""))
        BSLASH = "\"
        for (i=1; i<=strlen(to); i++) {
            if ((t=substr(to,i,1))==BSLASH) {
                i++
                if ((t=substr(to,i,1))==BSLASH) s = s + BSLASH      // "\\"
                else if ((j=strtoreal(t))<.)    s = s + sub[j+1]    // "\#"
                else                            s = s + BSLASH + t
            }
            else s = s + t
        }
        return(subinstr(s0, sub[1], s, 1))
    }
    else return(s0)
}

end


*! {smcl}
*! {marker mm_invtokens}{bf:mm_invtokens.mata}{asis}
*! version 1.0.2  Ben Jann  05jan2008
version 9.0
mata:

string scalar mm_invtokens(string vector In, | real scalar noclean)
{
    string scalar Out
    string scalar tok
    real scalar i

    if (args()<2) noclean = 0
    Out = ""
    for (i=1; i<=length(In); i++) {
        tok = In[i]
        if (noclean)                 tok = "`" + `"""' + tok + `"""' + "'"
        else if (strpos(tok, `"""')) tok = "`" + `"""' + tok + `"""' + "'"
        else if (strpos(tok, " "))   tok = `"""' + tok + `"""'
        else if (tok=="")            tok = `"""' + `"""'
        if (i>1) tok = " " + tok
        Out = Out + tok
    }
    return(Out)
}
end


*! {smcl}
*! {marker mm_realofstr}{bf:mm_realofstr.mata}{asis}
*! version 1.0.2, Ben Jann, 24mar2006
version 9.0
mata:

real matrix mm_realofstr(string matrix S)
{
	real matrix R
	string scalar tmp
	real scalar i, j

	R = J(rows(S), cols(S), .)
	tmp = st_tempname()
	for (i=1; i<=rows(S); i++) {
		for (j=1; j<=cols(S); j++) {
			stata("scalar " + tmp + "=real(`" + `"""' + S[i,j] + `"""' + "')")
			R[i,j] = st_numscalar(tmp)
		}
	}
	stata("capture scalar drop " + tmp)
	return(R)
}
end


*! {smcl}
*! {marker mm_strexpand}{bf:mm_strexpand.mata}{asis}
*! version 1.0.4, Ben Jann, 30apr2007
version 9.0
mata:

string scalar mm_strexpand(string scalar s, string vector slist,
 | string scalar def, real scalar unique, string scalar errtxt)
{
	real scalar   err
	string scalar res

	if (args()<5) errtxt = `"""' + s + `"" invalid"'
	if (args()<4) unique = 0
	err = _mm_strexpand(res, s, slist, def, unique)
	if (err) _error(err, errtxt)
	return(res)
}

real scalar _mm_strexpand(res, string scalar s,
 string vector slist, | string scalar def, real scalar unique)
{
	real scalar i, l, match

	if (s=="") {
		res = def
		return(0)
	}
	if (args()<5) unique = 0
	l = strlen(s)
	if (unique) {
		match = 0
		for (i=1; i<=length(slist); i++) {
			if (s==substr(slist[i], 1, l)) {
				if (match) return(3498)
				match = i
			}
		}
		if (match) {
			res = slist[match]
			return(0)
		}
	}
	else {
		for (i=1; i<=length(slist); i++) {
			if (s==substr(slist[i], 1, l)) {
				res = slist[i]
				return(0)
			}
		}
	}
	return(3499)
}
end


*! {smcl}
*! {marker mm_matlist}{bf:mm_matlist.mata}{asis}
*! version 1.0.4  31aug2007  Ben Jann
version 9.2
mata:

void mm_matlist(
    real matrix X,
    | string matrix fmt,
      real scalar b,
      string vector rstripe,
      string vector cstripe,
      string scalar rtitle,
      transmorphic matrix res
    )
{
    real scalar             i, j, rw, lw, j0, j1, l, r, save
    real rowvector          wd
    string rowvector        sfmt
    pointer(real scalar) scalar     fi, fj
    pointer(string vector) scalar   rstr, cstr

    if (args()<2) fmt = "%9.0g"
    if (args()<3) b = 1
    save = (args()==7)
    //wd = strlen(sprintf(fmt,-1/3))

    if ((rows(fmt)!=1 & rows(fmt)!=rows(X)) |
        (cols(fmt)!=1 & cols(fmt)!=cols(X))) _error(3200)
    if (!anyof((0,1,2,3),b)) _error(3300)

    fi = (rows(fmt)!=1 ? &i : &1)
    fj = (cols(fmt)!=1 ? &j : &1)

    wd = J(1,cols(X),0)
    for (i=1;i<=rows(X);i++) {
        for (j=1;j<=cols(X);j++) {
            wd[j] = max((wd[j],strlen(sprintf(fmt[*fi,*fj],X[i,j]))))
        }
    }

    if (length(X)==0) return

    if (length(rstripe)<1)  rstr = &strofreal(1::rows(X))
    else rstr = &rstripe
    if (length(cstripe)<1)  cstr = &strofreal(1..cols(X))
    else if (rows(cstripe)!=1)  cstr = &(cstripe')
    else                        cstr = &cstripe

    if ((*rstr!="" & length(*rstr)!=rows(X)) |
        (*cstr!="" & length(*cstr)!=cols(X))) _error(3200)

    if (*rstr!="")  rw = max((max(strlen(*rstr)),strlen(rtitle)))
    else            rw = -1
    if (*cstr!="")  wd = colmax((strlen(*cstr)\wd)) :+ 2
    else            wd = wd :+ 2

    sfmt = "%" :+ strofreal(wd) :+ "s"

    if (save) {
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, 1, cols(X), 1, 1, res)
        return
    }

    lw = st_numscalar("c(linesize)") - ((2+rw+2*(b>0))+2*(b==1))
    if ((sum(wd)+cols(X))<=lw) {
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, 1, cols(X), 1, 1)
        return
    }

    j0 = 1
    l = 1
    r = 0
    while (j0<=cols(X)) {
        j1 = j0
        while (sum(wd[|j0\j1|]:+1)<=lw) {
            if (++j1>cols(X)) {
                r = 1
                break
            }
        }
        j1 = max((j1-1,j0))
        _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle,
            wd, rw, sfmt, j0, j1, l, r)
        l = 0
        j0 = j1+1
    }
}

string colvector mm_smatlist(
    real matrix X,
    | string matrix fmt,
      real scalar b,
      string vector rstripe,
      string vector cstripe,
      string scalar rtitle
    )
{
    string colvector res

    if (args()<2) fmt = "%9.0g"
    if (args()<3) b = 1
    if (args()<4) rstripe = J(0,1,"")
    if (args()<5) cstripe = J(1,0,"")
    if (args()<5) rtitle = ""
    mm_matlist(X, fmt, b, rstripe, cstripe, rtitle, res)
    return(res)
}

void _mm_matlist(
    real matrix X,
    string matrix fmt,
    real scalar b,
    string vector rstripe,
    string vector cstripe,
    string scalar rtitle,
    real rowvector wd,
    real scalar rw,
    string rowvector sfmt,
    real scalar j0,
    real scalar j1,
    real scalar l,
    real scalar r,
    | transmorphic matrix res
    )
{
    real scalar          i, j, di, ri
    string scalar        rfmt, line
    pointer(real scalar) scalar fi, fj

    di = args()<14
    fi = (rows(fmt)!=1 ? &i : &1)
    fj = (cols(fmt)!=1 ? &j : &1)
    rfmt = "%"+strofreal(rw)+"s"

    if (di==0) {
        ri = 0
        res = J(rows(X)+(1+(b==3))*(cstripe!="")+(b>0)+(b==1)
            +(b==3)+(l==0&(b==0|b==2)),1,"")
    }

    if (l==0 & (b==0 | b==2)) {
        if (di) printf("\n")
        else res[++ri] = ""
    }
    if (cstripe!="") {
        if (b==3) {
            line = sprintf("{txt}" + "{hline " + strofreal(2+rw+1) +  "}{c TT}"
                + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)-1) + "}")
            if (di) display(line)
            else res[++ri] = line
        }
        line = sprintf("{txt}  ")
        if (rstripe!="") line = line + sprintf(rfmt+" ", rtitle)
        if (b>0) line = line + sprintf((b==1 ? " " : "{c |}"))
        for (j=j0;j<=j1;j++) {
            line = line + sprintf(sfmt[j], cstripe[j])
            if (j<j1) line = line + sprintf(" ")
        }
        if (di) display(line)
        else res[++ri] = line
    }
    if (b) {
        line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " +
            (l ? "{c TLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c +}")
            + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" +
            (b==1 ? (r ? "{c TRC}" : "") : ""))
        if (di) display(line)
        else res[++ri] = line
    }
    for (i=1;i<=rows(X);i++) {
        line = sprintf("{txt}  ")
        if (rstripe!="") line = line + sprintf(rfmt+" ", rstripe[i])
        if (b>0) line = line + sprintf((l | b!=1 ? "{c |}" : " "))
        line = line + sprintf("{res}")
        for (j=j0;j<=j1;j++) {
            line = line + sprintf(sfmt[j],sprintf(fmt[*fi,*fj], X[i,j]))
            if (j<j1) line = line + sprintf(" ")
        }
        if (b==1 & r) line = line + sprintf("  {txt}{c |}")
        if (di) display(line)
        else res[++ri] = line
    }
    if (b==1|b==3) {
        line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " +
            (l ? "{c BLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c BT}")
            + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" +
            (b==1 ? (r ? "{c BRC}" : "") : ""))
        if (di) display(line)
        else res[++ri] = line
    }
}

end


*! {smcl}
*! {marker mm_insheet}{bf:mm_insheet.mata}{asis}
*! version 2.0.1, Ben Jann, 18mar2006
*! based on cat(), version 2.0.0  23jan2006
version 9.0
mata:

string matrix mm_insheet(string scalar filename,
 | string scalar del, real scalar uline1, real scalar line2)
{
	real scalar    fh, fpos, line1, i, r, c, j, delpos, ji
	string matrix  EOF
	string matrix  res
	string scalar  line

	/* ------------------------------------------------------------ */
	/* setup */
	if (del=="") del = char(9)
	fh  = fopen(filename, "r")
	EOF = J(0, 0, "")
	line1 = floor(uline1)
	if (line1<1 | line1>=.) line1 = 1

	/* ------------------------------------------------------------ */
	/* nothing to read case  */
	if (line1>line2) {
		fclose(fh)
		return(J(0, 0, ""))
	}

	/* ------------------------------------------------------------ */
	 /* skip opening lines  */
	for (i=1; i<line1; i++) {
		if (fget(fh)==EOF) {
			fclose(fh)
			return(J(0, 0, ""))
		}
	}
	fpos = ftell(fh)

	/* ------------------------------------------------------------ */
	/* count lines and columns  */
	j = 0
	for (i=line1; i<=line2; i++) {
		if ((line=fget(fh))==EOF) break
		ji = 1
		while (delpos = strpos(line, del)) {
			ji++
			line = substr(line, delpos+1, .)
		}
		if (ji>j) j = ji
	}
	res = J(r = i-line1, c = j, "")

	/* ------------------------------------------------------------ */
	/* read lines */
	fseek(fh, fpos, -1)
	for (i=1; i<=r; i++) {
		if ((line=fget(fh))==EOF) {
			/* unexpected EOF -- file must have changed */
			fclose(fh)
			if (--i) return(res[|1,1 \ i,c|])
			return(J(0,c,""))
		}
		for (j=1; j<=c; j++) {
			delpos = strpos(line, del)
			if (delpos==0) {
				res[i,j] = line
				break
			}
			res[i,j] = substr(line, 1, delpos-1)
			line = substr(line, delpos+1, .)
		}
	}
	fclose(fh)
	return(res)
}

end


*! {smcl}
*! {marker mm_infile}{bf:mm_infile.mata}{asis}
*! version 2.0.2, Ben Jann, 18mar2006
*! based on cat(), version 2.0.0  23jan2006
version 9.0
mata:

string matrix mm_infile(string scalar filename,
 | real scalar uline1, real scalar line2)
{
	real scalar    fh, fpos, line1, i, r, c, j
	string matrix  EOF
	string matrix  res
	string scalar  line

	/* ------------------------------------------------------------ */
	/* setup */
	fh  = fopen(filename, "r")
	EOF = J(0, 0, "")
	line1 = floor(uline1)
	if (line1<1 | line1>=.) line1 = 1

	/* ------------------------------------------------------------ */
	/* nothing to read case  */
	if (line1>line2) {
		fclose(fh)
		return(J(0, 0, ""))
	}

	/* ------------------------------------------------------------ */
	 /* skip opening lines  */
	for (i=1; i<line1; i++) {
		if (fget(fh)==EOF) {
			fclose(fh)
			return(J(0, 0, ""))
		}
	}
	fpos = ftell(fh)

	/* ------------------------------------------------------------ */
	/* count lines and columns  */
	j = 0
	for (i=line1; i<=line2; i++) {
		if ((line=fget(fh))==EOF) break
		line = tokens(line)
		if (cols(line)>j) j = cols(line)
	}
	res = J(r = i-line1, c = j, "")

	/* ------------------------------------------------------------ */
	/* read lines */
	fseek(fh, fpos, -1)
	for (i=1; i<=r; i++) {
		if ((line=fget(fh))==EOF) {
			/* unexpected EOF -- file must have changed */
			fclose(fh)
			if (--i) return(res[|1,1 \ i,c|])
			return(J(0,c,""))
		}
		line = tokens(line)
		res[|i,1 \ i,cols(line)|] = line
	}
	fclose(fh)
	return(res)
}

end


*! {smcl}
*! {marker mm_outsheet}{bf:mm_outsheet.mata}{asis}
*! version 1.0.4, Ben Jann, 14apr2006
version 9.0
mata:

void mm_outsheet(string scalar fn, string matrix s,
 | string scalar mode0, string scalar del)
{
	string scalar line, mode, m
	real scalar i, j, fh

	if (args()<4) del = char(9)
	mode = mm_strexpand(mode0, ("append", "replace"))
	m = "w"
	if (mode=="replace") unlink(fn)
	else if (mode=="append") m = "a"
	fh = fopen(fn, m)
	for (i=1; i<=rows(s); i++) {
		line = J(1,1,"")
		for (j=1; j<=cols(s); j++) {
			line = line + s[i,j]
			if (j<cols(s)) line = line + del
		}
		fput(fh, line)
	}
	fclose(fh)
}
end


*! {smcl}
*! {marker mm_callf}{bf:mm_callf.mata}{asis}
*! version 1.0.1, Ben Jann, 16jun2006
version 9.2

local o1 "o1"
local po1 "*p.o1"
forv i=2/10 {
	local o`i' "`o`=`i'-1'', o`i'"
	local po`i' "`po`=`i'-1'', *p.o`i'"
}

mata:

struct mm_callf_o10 {
	pointer(function) scalar f
	real scalar              n
	pointer scalar           `o10'
}

struct mm_callf_o10 scalar mm_callf_setup(
 pointer(function) scalar f, real scalar n, | `o10')
{
	struct mm_callf_o10 scalar p

	p.f  = f
	p.n  = n
	p.o1 = &o1 ; p.o2 = &o2 ; p.o3 = &o3 ; p.o4 = &o4 ; p.o5 = &o5
	p.o6 = &o6 ; p.o7 = &o7 ; p.o8 = &o8 ; p.o9 = &o9 ; p.o10 = &o10
	return(p)
}

transmorphic mm_callf(struct mm_callf_o10 p, | `o5')
{
	if (args()==1) return(_mm_callf0(p))
	if (args()==2) return(_mm_callf1(p, `o1'))
	if (args()==3) return(_mm_callf2(p, `o2'))
	if (args()==4) return(_mm_callf3(p, `o3'))
	if (args()==5) return(_mm_callf4(p, `o4'))
	               return(_mm_callf5(p, `o5'))
}

transmorphic _mm_callf0(struct mm_callf_o10 p)
{
	if (p.n<=0) return( (*p.f)() )
	if (p.n==1) return( (*p.f)(`po1') )
	if (p.n==2) return( (*p.f)(`po2') )
	if (p.n==3) return( (*p.f)(`po3') )
	if (p.n==4) return( (*p.f)(`po4') )
	if (p.n==5) return( (*p.f)(`po5') )
	if (p.n==6) return( (*p.f)(`po6') )
	if (p.n==7) return( (*p.f)(`po7') )
	if (p.n==8) return( (*p.f)(`po8') )
	if (p.n==9) return( (*p.f)(`po9') )
	            return( (*p.f)(`po10') )
}

transmorphic _mm_callf1(struct mm_callf_o10 p, `o1')
{
	if (p.n<=0) return( (*p.f)(`o1') )
	if (p.n==1) return( (*p.f)(`o1', `po1') )
	if (p.n==2) return( (*p.f)(`o1', `po2') )
	if (p.n==3) return( (*p.f)(`o1', `po3') )
	if (p.n==4) return( (*p.f)(`o1', `po4') )
	if (p.n==5) return( (*p.f)(`o1', `po5') )
	if (p.n==6) return( (*p.f)(`o1', `po6') )
	if (p.n==7) return( (*p.f)(`o1', `po7') )
	if (p.n==8) return( (*p.f)(`o1', `po8') )
	if (p.n==9) return( (*p.f)(`o1', `po9') )
	            return( (*p.f)(`o1', `po10') )
}

transmorphic _mm_callf2(struct mm_callf_o10 p, `o2')
{
	if (p.n<=0) return( (*p.f)(`o2') )
	if (p.n==1) return( (*p.f)(`o2', `po1') )
	if (p.n==2) return( (*p.f)(`o2', `po2') )
	if (p.n==3) return( (*p.f)(`o2', `po3') )
	if (p.n==4) return( (*p.f)(`o2', `po4') )
	if (p.n==5) return( (*p.f)(`o2', `po5') )
	if (p.n==6) return( (*p.f)(`o2', `po6') )
	if (p.n==7) return( (*p.f)(`o2', `po7') )
	if (p.n==8) return( (*p.f)(`o2', `po8') )
	if (p.n==9) return( (*p.f)(`o2', `po9') )
	            return( (*p.f)(`o2', `po10') )
}

transmorphic _mm_callf3(struct mm_callf_o10 p, `o3')
{
	if (p.n<=0) return( (*p.f)(`o3') )
	if (p.n==1) return( (*p.f)(`o3', `po1') )
	if (p.n==2) return( (*p.f)(`o3', `po2') )
	if (p.n==3) return( (*p.f)(`o3', `po3') )
	if (p.n==4) return( (*p.f)(`o3', `po4') )
	if (p.n==5) return( (*p.f)(`o3', `po5') )
	if (p.n==6) return( (*p.f)(`o3', `po6') )
	if (p.n==7) return( (*p.f)(`o3', `po7') )
	if (p.n==8) return( (*p.f)(`o3', `po8') )
	if (p.n==9) return( (*p.f)(`o3', `po9') )
	            return( (*p.f)(`o3', `po10') )
}

transmorphic _mm_callf4(struct mm_callf_o10 p, `o4')
{
	if (p.n<=0) return( (*p.f)(`o4') )
	if (p.n==1) return( (*p.f)(`o4', `po1') )
	if (p.n==2) return( (*p.f)(`o4', `po2') )
	if (p.n==3) return( (*p.f)(`o4', `po3') )
	if (p.n==4) return( (*p.f)(`o4', `po4') )
	if (p.n==5) return( (*p.f)(`o4', `po5') )
	if (p.n==6) return( (*p.f)(`o4', `po6') )
	if (p.n==7) return( (*p.f)(`o4', `po7') )
	if (p.n==8) return( (*p.f)(`o4', `po8') )
	if (p.n==9) return( (*p.f)(`o4', `po9') )
	            return( (*p.f)(`o4', `po10') )
}

transmorphic _mm_callf5(struct mm_callf_o10 p, `o5')
{
	if (p.n<=0) return( (*p.f)(`o5') )
	if (p.n==1) return( (*p.f)(`o5', `po1') )
	if (p.n==2) return( (*p.f)(`o5', `po2') )
	if (p.n==3) return( (*p.f)(`o5', `po3') )
	if (p.n==4) return( (*p.f)(`o5', `po4') )
	if (p.n==5) return( (*p.f)(`o5', `po5') )
	if (p.n==6) return( (*p.f)(`o5', `po6') )
	if (p.n==7) return( (*p.f)(`o5', `po7') )
	if (p.n==8) return( (*p.f)(`o5', `po8') )
	if (p.n==9) return( (*p.f)(`o5', `po9') )
	            return( (*p.f)(`o5', `po10') )
}

end


*! {smcl}
*! {marker mm_version}{bf:mm_version.mata}{asis}
*! version 2.0.0  27jan2022  Ben Jann

version 9.1
mata:

real scalar mm_version() return(200)

end
